**Engineering Applications of Computational Methods 13**

Chih-Yung Wen Yazhong Jiang Lisong Shi

# Space–Time Conservation Element and Solution Element Method

Advances and Applications in Engineering Sciences

## **Engineering Applications of Computational Methods**

Volume 13

#### **Series Editors**

Liang Gao, State Key Laboratory of Digital Manufacturing Equipment and Technology, Huazhong University of Science and Technology, Wuhan, Hubei, China

Akhil Garg, School of Mechanical Science and Engineering, Huazhong University of Science and Technology, Wuhan, Hubei, China

The book series Engineering Applications of Computational Methods addresses the numerous applications of mathematical theory and latest computational or numerical methods in various fields of engineering. It emphasizes the practical application of these methods, with possible aspects in programming. New and developing computational methods using big data, machine learning and AI are discussed in this book series, and could be applied to engineering fields, such as manufacturing, industrial engineering, control engineering, civil engineering, energy engineering and material engineering.

The book series Engineering Applications of Computational Methods aims to introduce important computational methods adopted in different engineering projects to researchers and engineers. The individual book volumes in the series are thematic. The goal of each volume is to give readers a comprehensive overview of how the computational methods in a certain engineering area can be used. As a collection, the series provides valuable resources to a wide audience in academia, the engineering research community, industry and anyone else who are looking to expand their knowledge of computational methods.

This book series is indexed in both the **Scopus** and **Compendex** databases.

## Space–Time Conservation Element and Solution Element Method

Advances and Applications in Engineering Sciences

Chih-Yung Wen Department of Aeronautical and Aviation Engineering, Faculty of Engineering The Hong Kong Polytechnic University Hong Kong, China

Lisong Shi Department of Aeronautical and Aviation Engineering, Faculty of Engineering The Hong Kong Polytechnic University Hong Kong, China

Yazhong Jiang Hubei Key Laboratory of Theory and Application of Advanced Materials Mechanics, School of Science Wuhan University of Technology Wuhan, Hubei, China

ISSN 2662-3366 ISSN 2662-3374 (electronic) Engineering Applications of Computational Methods ISBN 978-981-99-0875-2 ISBN 978-981-99-0876-9 (eBook) https://doi.org/10.1007/978-981-99-0876-9

© The Editor(s) (if applicable) and The Author(s) 2023. This book is an open access publication.

**Open Access** This book is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this book are included in the book's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the book's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

The use of general descriptive names, registered names, trademarks, service marks, etc. in this publication does not imply, even in the absence of a specific statement, that such names are exempt from the relevant protective laws and regulations and therefore free for general use.

The publisher, the authors, and the editors are safe to assume that the advice and information in this book are believed to be true and accurate at the date of publication. Neither the publisher nor the authors or the editors give a warranty, expressed or implied, with respect to the material contained herein or for any errors or omissions that may have been made. The publisher remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

This Springer imprint is published by the registered company Springer Nature Singapore Pte Ltd. The registered company address is: 152 Beach Road, #21-01/04 Gateway East, Singapore 189721, Singapore

### **Preface**

Methods for numerical simulations in engineering sciences keep updating every day, and naturally the state-of-the-art advancement may not be incorporated into the well-known textbooks. This short book introduces the Conservation Element and Solution Element (CESE) method, which has a relatively short history. The CESE method provides time-accurate simulations of physical systems and overcomes the difficulties in capturing discontinuities, by combining many innovative considerations in a nontraditional manner. This book elaborates on how the CESE method works and presents a comprehensive survey of the recent CESE developments and applications.

We would like to appreciate Dr. Ben Guan at Harbin Engineering University and Dr. Zijian Zhang at The Hong Kong Polytechnic University for their direct contributions to Chaps. 7 and 8, respectively. We are grateful to Prof. Hua Shen at University of Electronic Science and Technology of China for his important contribution to the CESE method, and also for his unselfish assistance to us in preparing this book. All group members in High-Speed Thermo-Fluid and MAV/UAV Laboratory at The Hong Kong Polytechnic University are deeply appreciated for their efforts and supports. We would like to thank Prof. De-Liang Zhang at Chinese Academy of Sciences, who has provided us invaluable inspiration to improve and apply the CESE method. We are also appreciative to Prof. Vigor Yang at Georgia Institute of Technology, USA, for his encouragement to us to write this book on the CESE method. We would like to acknowledge the support from Research Institute for Sports Science and Technology, The Hong Kong Polytechnic University (project code: 1-CD6N).

Hong Kong, China Wuhan, China Hong Kong, China Chih-Yung Wen Yazhong Jiang Lisong Shi

## **Contents**




## **Chapter 1 Introduction**

With the rapid development of electronic computers, numerical computation has become an important paradigm of scientific discovery as well as a powerful tool for engineering research. Solving complex problems in a computational fashion is more than applying theories, equations, and formulas. Computational methods, also called algorithms or schemes, have strong influences on the outcomes of computations.

The space–time conservation element and solution element (CESE) method is aimed to numerically solve equations of conservation laws in various physical systems. The major concern herein will be the computational fluid dynamics (CFD), which is a hot area in scientific and engineering researches. Nevertheless, it may help at the outset to recognize that the conservation laws in CFD problems have a lot in common with conservation laws in other fields, such as acoustics, solid dynamics, electromagnetics, and magnetohydrodynamics. The physics may be different, but the mathematics are similar. For instance, physics involving dynamical evolution of waves and discontinuities are usually modelled by time-dependent nonlinear hyperbolic partial differential equations (PDEs). Some CFD problems happen to be representative of such physical problems in a wider context.

#### **1.1 Background**

Most of the key achievements in conventional CFD have been incorporated into the computational procedure of the finite volume method (FVM). Indeed, the FVM is well established and widely used, and therefore it is sometimes regarded as a mature technique. There are two important steps in the FVM, namely the reconstruction and the evolution [1, 2]. In the reconstruction step, one needs to locally approximate the flow field with simple functions such as polynomials, and then interpolate the values of flow quantities at the cell interfaces by utilizing the cell-average values. In the evolution step, the numerical flux at each cell interface is determined by some kind

<sup>©</sup> The Author(s) 2023

C.-Y. Wen et al., *Space–Time Conservation Element and Solution Element Method*, Engineering Applications of Computational Methods 13, https://doi.org/10.1007/978-981-99-0876-9\_1

of flux technique, and then the time integration of the conservation equations can be performed.

A major category of CFD problems face the challenges brought by (1) the coexistence of smooth regions and discontinuities in the flow (e.g., shock waves and material interfaces) and (2) the wide spectrum of wave numbers/frequencies in the flow (e.g., acoustic waves and turbulence). Due to these physical phenomena, the CFD solver is required to well handle the numerical dissipation, which leads to an endless debate in the CFD community. First, for the robustness in capturing shock waves, there must be a certain amount of numerical dissipation to suppress the spurious oscillations. However, a low numerical dissipation is very favourable for the accuracy and resolution of the viscous layers, contact surfaces, and small-scale structures in the flow, which would be smeared out by excessive dissipation.

To achieve the low but controllable numerical dissipation, the research works on FVM are mostly dedicated to two approaches: (1) developing high-order reconstruction techniques in conjunction with nonlinear limiters and (2) improving the flux solvers (e.g., approximate Riemann solvers). In both ways, numerous important advances have been obtained.

In spite of these successes in the FVM development, several shortcomings may arise. (1) A conventional FVM implementation highly relies on the selections of special techniques, such as the limiter and the Riemann solver. Actually, each special technique has its own pros and cons. There are many alternatives for each option, but none of them can be proven to be optimal and general. (2) The coupling of space and time is usually overlooked and may not be guaranteed in the numerical solution. (3) The multi-dimensionality can be questionable because most flux solvers are only based on one-dimensional physics. (4) The popular high-order FVM schemes are usually not compact, i.e., large stencils are used. Under this background, the CESE method was proposed to overcome the difficulties in the conventional FVM to some extent, and it became a new member in the family of FVM.

#### **1.2 History of the CESE Method**

The space–time conservation element and solution element method was initiated by Sin-Chung Chang of NASA Glenn Research Centre and his collaborators for aerodynamic computation in early 1990s [3, 4]. There were several motivations to propose the CESE method. First, they wanted to construct multidimensional and space–time unified CFD schemes, without using dimensional splitting or separated treatments of space discretization and time discretization. Second, they believed a CESE solver should be built from a non-dissipative core scheme so that numerical dissipation can be controlled effectively, dynamically, and even actively. Third, they attempted to avoid using the Riemann solver in the CESE method.

Around 1994, the CESE project was approved and supported by NASA Glenn Research Centre. The work of the research group showed that this time-accurate scheme possesses low numerical dissipation, which is valuable in CFD simulation. For the first time, the transonic resonance in a convergent-divergent nozzle with frequency-staging effect was successfully simulated by the CESE scheme [5, 6]. In this case, the blind prediction by the CESE scheme matched the experimental data. Therefore, the CESE method is proved to be capable of simultaneously simulating multidimensional unsteady shock waves and acoustic waves with high accuracy.

Owing to its accuracy and robustness, the CESE method is employed in the CFD module of the simulation software LS-DYNA [7]. The Jacobs Technology Inc. (designer of the NASA's hypersonic flow test facilities) developed its own in-house CESE code called JUSTUS (Jacobs Unified Space–Time Unstructured Solver) for hypersonic flow simulations. Another in-house CESE code, called "ez4d" software, was developed by the NASA Langley Research Centre. This is a time-accurate 3D Navier–Stokes flow solver on unstructured meshes [8, 8].

#### **1.3 Main Features of the CESE Method**

The CESE method is a special finite-volume-type numerical method, with the aim of solving the governing equations of fluid dynamics and other conservation laws in various physical systems. The CESE method possesses many features such as the unified treatment for space and time, the fully discrete one-step explicit scheme, and the highly compact stencil. Without enlarging the stencil or adding stages of time integration, the CESE method can achieve arbitrary high-order accuracy for both space and time.

The essential ingredients in the CESE method include: (1) the spatial derivatives of physical variables are stored as independent unknowns, in addition to the physical variables themselves. In every time step, these derivatives are updated by a specially designed procedure. (2) a staggered mesh and a staggered time-marching strategy are employed. (3) the interior structure within each solution element is built with the Taylor expansion. (4) the time-marching approach is based on the Cauchy–Kowalewski procedure.

Through a combination of the above techniques, the CESE method has low dissipation and high compactness. When applied to the simulations of complex physical processes, the CESE method can catch shock waves, contact discontinuities, fine structures, and small disturbances, with high resolution and strong robustness. Therefore, the CESE method demonstrates good performances in the numerical simulations of wave-propagation problems (e.g., shock waves, acoustic waves, detonation waves, stress waves in solid, and the electromagnetic waves), interfacial instabilities, as well as the interactions of gaseous and liquid phases. In many research areas including high-speed aerodynamics, shock dynamics, detonation, aeroacoustics, and solid dynamics, the CESE method proves to be suitable and shows a good development prospect.

#### **1.4 Outline of the Book**

The remainder of this book is organized as follows. First of all, the non-dissipative core scheme of the CESE method is introduced in Chap. 2, with detailed descriptions of the basic concepts. Then, the practical shock-capturing CESE schemes with numerical dissipation including the classical *a*–α scheme, the Courant number insensitive scheme, and the recently proposed upwind CESE schemes are presented in Chap. 3. Furthermore, Chaps. 4 and 5 extend the CESE method to multidimensional and high-order versions. In Chap. 6, numerical properties of various CESE schemes are analysed, along with comparisons to other numerical schemes. Chapters 7–9 provide an overview of the applications of the CESE method, most of which are quite relevant to the aerospace engineering. Finally, a summary and an outlook of the CESE method are given in Chap. 10.

#### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **Chapter 2 Non-dissipative Core Scheme of CESE Method**

This chapter is devoted to demonstrating the basic ideas in the CESE method. These ideas include the adoption of a space–time integral form of governing equations as the starting point of scheme construction, as well as the introduction of conservation element (CE) and solution element (SE) in the discretization of space–time domain. Then, the non-dissipative core scheme of the CESE method will be presented in detail.

#### **2.1 Space–Time Integral Form of Governing Equations**

Governing equations for a specific physical problem can be expressed in both differential and integral forms. Usually, a differential equation or a set of differential equations can clearly describe the evolution of a physical system. The differential forms of governing equations are prevalent in textbooks and theoretical research works, due to their compactness in writing and the maturity in the mathematical analysis on differential equations.

In the area of numerical simulation, the finite difference method is directly applied to the differential equations, while the finite volume method and finite element method require the governing equations to be cast into integral forms before numerical treatments. In fact, the performance of a numerical method will be influenced by the form of equations that is used in the numerical method, even though different forms of governing equations can be mathematically equivalent. As will be shown in this book, a major feature of the CESE method is the adoption of a space–time integral form of the governing equations, in which time and space are treated on the same footing.

For illustrative purposes, we consider the compressible Euler equations for 2D planar flows, which can be written in a differential form as

8 2 Non-dissipative Core Scheme of CESE Method

$$
\frac{
\partial \mathbf{U}
}{
\partial t
} + \frac{
\partial \mathbf{F}
}{
\partial x
} + \frac{
\partial \mathbf{G}
}{
\partial y
} = \mathbf{0},
\tag{2.1}
$$

where *U* denotes the vector of conserved variables such that

$$\mathbf{U} = \begin{pmatrix} \rho \\ \rho u \\ \rho v \\ E \end{pmatrix},\tag{2.2}$$

where *F* and *G* are the inviscid fluxes of the form

$$\mathbf{F} = \begin{pmatrix} \rho u \\ \rho u^2 + p \\ \rho uv \\ (E+p)u \end{pmatrix}, \mathbf{G} = \begin{pmatrix} \rho v \\ \rho uv \\ \rho v^2 + p \\ (E+p)v \end{pmatrix} \tag{2.3}$$

In Eqs. (2.2) and (2.3), ρ is the density of the fluid, *u* and *v* are the *x*-component and *y*-component of the flow velocity, respectively, *p* is the static pressure, and *E* is the total energy per unit volume of the fluid. Note that once the equation of state (EOS) is provided, Eq. (2.1) becomes a closed set of equations. Here, flux vectors *F* and *G* are regarded as functions of conserved vector *U*. A set of equations in the form of Eq. (2.1) is called a system of conservation laws.

By using the gradient operator in 2D physical space

$$\nabla \equiv \left(\frac{\partial}{\partial \mathbf{x}}, \frac{\partial}{\partial \mathbf{y}}\right), \tag{2.4}$$

Equation (2.1) can also be written in a divergence form

$$
\frac{
\partial \mathbf{U}
}{
\partial t
} + \nabla \cdot \mathbf{H} = \mathbf{0},
\tag{2.5}
$$

where *H* is a matrix composed of spatial flux vectors *F* and *G*:

$$\mathbf{H} = (\mathbf{F}, \mathbf{G})\tag{2.6}$$

Therefore, the divergence form of Euler equation states that at time *t* and point (*x*, *y*), the temporal rate of change of conserved quantities plus the divergence of the spatial flux of these quantities must be zero.

Consider a control volume denoted by *V* in 2D physical space, as shown in Fig. 2.1, with the surface *S*(*V*) and the unit outward normal vector *n* on the surface. Integrating Eq. (2.5) over the control volume *V* and applying the Gauss's theorem, one obtains

**Fig. 2.1** Schematic of a control volume in two-dimensional physical space

$$\frac{\partial}{\partial t} \left( \int\_{V} \mathbf{U} \, \mathrm{d}V \right) + \oint\_{S(V)} \mathbf{H} \cdot \mathbf{n} \, \mathrm{d}S = \mathbf{0},\tag{2.7}$$

which can be interpreted as follows: the temporal rate of change of conserved quantities within the volume *V* must equal to the net inward flow rate of these quantities through the surface *S*(*V*). This integral form of Euler Eq. (2.7) is the starting point of the well-known finite volume method (FVM) in computational fluid dynamics.

With a unified treatment of time and physical space, we can define *x*, *y,* and *t* as coordinates in a three-dimensional Euclidean space, called the space–time domain. On this basis, the definition of gradient operator can be extended as

$$\nabla \equiv \left(\frac{\partial}{\partial x}, \frac{\partial}{\partial y}, \frac{\partial}{\partial t}\right) \tag{2.8}$$

and Eq. (2.1) can be written in a divergence-free form as

$$\nabla \cdot \mathbf{h} = \mathbf{0},\tag{2.9}$$

where the matrix *h* is composed of both flux vectors and the conserved vector *U*:

$$\mathbf{h} = (\mathbf{F}, \mathbf{G}, \mathbf{U}) \tag{2.10}$$

By applying Gauss's theorem to the integration of Eq. (2.9) over an arbitrary control volume *V* in the 3D space–time domain as shown in Fig. 2.2, the Euler equation is eventually formulated as

$$\oint\_{S(V)} \mathbf{h} \cdot \mathbf{n} \, dS = \mathbf{0} \tag{2.11}$$

which means the balancing of flux is enforced for a space–time control volume. Clearly, this equation gives a direct description of the space–time conservation of mass, momentum, and energy in fluid flows, which faithfully preserves the original physical laws.

Equation (2.11) provides an example of the space–time integral form of governing equations for physical problems. It is worth noting that the space–time integral form and the traditional integral form like Eq. (2.7) state the same physics and must be equivalent to each other in mathematics, but they can lead to different numerical schemes. In contrast to the finite volume method, the CESE method numerically implements the space–time integral form of governing equations, emphasizing on the conservative property in the unity of time and physical space.

#### **2.2 Definitions of Conservation Elements and Solution Elements**

The CESE method starts with discretizing the space–time domain which is relevant to the computation. The discretization procedure generates the mesh for the physical space and determines the time-step size for the time-marching algorithm, similar to most CFD methods. The features of the CESE discretization include the construction of control volumes for the space–time integral form of the governing equation, the arrangement of the solution points, and the selection of the unknown variables to be calculated and stored at each solution point. All these features can be demonstrated by introducing two special concepts: the conservation element (CE) and the solution element (SE), which coin the algorithm accordingly.

To explain the CESE method in a simple way, we begin with the application of the CESE method to a 1D scalar conservation law of the form, without loss of generality

$$\frac{\partial u}{\partial t} + \frac{\partial f(u)}{\partial x} = 0.\tag{2.12}$$

With a uniform division of the 1D physical space and a constant time step, the space–time mesh is constructed as a 2D *x*–*t* plane as shown in Fig. 2.3 (solid lines). The spatial coordinate of the *j*-th mesh node is denoted by *xj*, with the corresponding cell centre position at *xj*+1/2 = (*xj* + *xj*+1)/2, and the cell size of Δ*x* = *xj*+1 − *xj*. In this

**Fig. 2.3** Mesh and arrangement of solution points for 1D CESE scheme

time-marching algorithm, each step between *tn*−1 and *tn* consists of two half steps: [*tn*−1, *tn*−1/2] and [*tn*−1/2, *tn*], with the time-step size of Δ*t* = *tn* − *tn*−1.

The unknown function *u*(*x*, *t*) is represented by the discrete values of *u* at a set of specific space–time points, called solution points which are shown as mesh nodes in Fig. 2.3—black circles for integer time levels {*t*0, *t*1, …, *tn*, …} and red squares (cell centres) for half-integer time levels {*t*1/2, *t*3/2, …, *tn*+1/2,…}. In other words, a staggered mesh is used at each intermediate time level. In the CESE scheme, both the unknown variable *u* and its spatial derivative *ux* will be calculated and stored at each solution point (e.g. point (*j*, *n*)):

$$\left(u\_j'' \equiv u(\mathbf{x}\_j, t\_n), \left(u\_x\right)\_j'' \equiv \frac{\partial u}{\partial x}(\mathbf{x}\_j, t\_n)\right) \tag{2.13}$$

The paths of information flow in a single CESE time step is illustrated in the schematic diagram in Fig. 2.4. As seen, a highly compact stencil in space and time is used in the CESE scheme. If a half time step it is treated as the basic iteration, the halfwidth of the symmetric stencil is then Δ*x*/2 and unknowns at point (*j*, *n*) only depend on the data stored at (*j* − 1/2, *n* − 1/2) and (*j* + 1/2, *n* − 1/2).

Based on the space–time mesh shown in Fig. 2.3, a set of small space–time elements, named conservation elements (CE), can be constructed, as demonstrated in Fig. 2.5. A CE is assigned to each solution point. For example, (CE) *n <sup>j</sup>*is denoted as the CE for point ( *j*, *n*), which is the rectangle with four vertices at the points (*j* − 1/2, *n* − 1/2), ( *j* + 1/2, *n* − 1/2), (*j* + 1/2, *n*), and (*j* − 1/2, *n*). Apparently, the CEs cover the whole space–time domain without overlap. It is noteworthy that

**Fig. 2.4** The information flow in one CESE time step with time-marching variables

the arrangement of CEs is staggered in two successive half time steps. Accordingly, the space–time integral form of Eq. (2.12) is numerically implemented within each conservation element and discrete equations for unknowns are established.

When performing the space–time integration of Eq. (2.12) over a CE, an important question arises: how to evaluate *u* and *f* along the boundary of the CE. This leads to the

**Fig. 2.5** Schematic of conservation elements (CEs) arrangement

**Fig. 2.6** Schematic of the solution element (SE) associated with solution point (*j*, *n*)

introduction of a solution element (SE) for each solution point, as shown in Fig. 2.6. (SE) *n <sup>j</sup>*composed of two line segments bisecting each other at point (*j*, *n*), forming a cross with four endpoints (*j* + 1/2, *n*), (*j*, *n* + 1/2), (*j* − 1/2, *n*), and (*j*, *n* − 1/2). For half-integer points, SEs can be defined in the same way. Note that conservation element (CE) *n <sup>j</sup>*is bounded by line segments belonging to three solution elements: (SE) *n <sup>j</sup>*, (SE) *n*−1/2 *<sup>j</sup>*−1/<sup>2</sup>, and (SE) *n*−1/2 *<sup>j</sup>*+1/<sup>2</sup>. In the CESE scheme, the functions *u*(*x*, *t*) and *f* (*x*, *t*) are assumed to be linear along each line segment of a SE and are approximated by first-order Taylor expansions about the centre of the SE. Specifically, in the solution element (SE) *n <sup>j</sup>*, *u* and *f* are constructed as

$$u(\mathbf{x},t) = u\_j^n + (u\_x)\_j^n(\mathbf{x} - \mathbf{x}\_j) + (u\_t)\_j^n(t - t\_n), \ (\mathbf{x}, t) \in (\mathbf{SE})\_j^n \tag{2.14}$$

$$f(\mathbf{x}, t) = f\_j^{\pi} + (f\_x)\_j^{\pi} (\mathbf{x} - \mathbf{x}\_j) + (f\_t)\_j^{\pi} (t - t\_n), \ (\mathbf{x}, t) \in (\mathbf{SE})\_j^{\pi} \tag{2.15}$$

where *ut* and *f t* are the temporal derivatives of *u* and *f* , respectively.

#### **2.3 Non-dissipative Core Scheme: a Scheme**

In this section, we present a non-dissipative CESE scheme, named the *a* scheme, to solve the 1D scalar conservation law, Eq. (2.12). A space–time flux vector is defined as

$$\mathbf{h} = (f, u) \tag{2.16}$$

where *f* and *u* can be regarded as the components of flux vector *h* in *x*-direction and *t*-direction, respectively. According to this definition, Eq. (2.12) can be converted into the space–time integral form of Eq. (2.11) by following the procedure in Sect. 2.1.

Consider a half step marching from time level *n*−1/2 to time level *n*. At the solution point (*j*, *n*), two independent unknowns *u<sup>n</sup> <sup>j</sup>*and (*ux* ) *n j* , need to be calculated simultaneously. Therefore, two algebraic equations need to be formulated by discretizing

the integral conservation law. For this purpose, the conservation element (CE) *n <sup>j</sup>*is split into two sub-elements: (CE−) *n <sup>j</sup>*and (CE+) *n <sup>j</sup>*. As shown in Fig. 2.7, (CE−) *n <sup>j</sup>*and (CE+) *n <sup>j</sup>*are the rectangle ACDE and the rectangle CBFD, respectively. Each edge of these rectangles belongs to one of the three SEs associated with (CE) *n j* .

Next, the space–time integral conservation law is implemented over each of the sub-CEs. Let the control volume *V* in Eq. (2.11) be (CE−) *n <sup>j</sup>*and (CE+) *n <sup>j</sup>*in turn, one can obtain

$$\oint\_{S(CE^{-}); \,} \mathbf{h} \cdot \mathbf{n} \,\mathrm{d}S = 0 \tag{2.17}$$

and

$$\oint\_{S(\mathbf{C}\mathbf{E}^+)\_{\neq}^\*} \mathbf{h} \cdot \mathbf{n} \,\mathrm{d}S = 0,\tag{2.18}$$

where *h* is the space–time flux vector expressed in Eq. (2.16) and *n* is the unit outward normal vector on the boundary of (CE−) *n <sup>j</sup>*or (CE+) *n j* .

As depicted in Fig. 2.8, the average values of *u*(*x*, *t*) on line segments DE, DF, AC, and BC are denoted by *UL*  \*, *UR*  \*, *UL*, and *UR*, respectively, while the average values of *f* (*x*, *t*) on line segments AE, BF, and CD are presented by *FL*, *FR*, and *FC*, respectively. Notably, CD is the interface between two sub-CEs. With this notation, integral Eqs. (2.17) and (2.18) can be expressed as

$$U\_L^\* \frac{\Delta \chi}{2} = U\_L \frac{\Delta \chi}{2} + (F\_L - F\_C) \frac{\Delta t}{2},\tag{2.19}$$

$$U\_R^\* \frac{\Delta \chi}{2} = U\_R \frac{\Delta \chi}{2} + (F\_C - F\_R) \frac{\Delta t}{2},\tag{2.20}$$

which explicitly state the balance of space–time flux in each sub-CE.

**Fig. 2.8** Space–time flux through the each boundary of sub-CEs

To proceed with the scheme construction, the time-marching variables (i.e. *u* and *ux*) at the solution points will be linked with *UL*  \*, *UR*  \*, *UL*, *UR*, *FL*, *FR*, and *FC*  in Eqs. (2.19) and (2.20), using the concept of the solution element. Recall that *u*  and *f* can be calculated by first-order Taylor expansions about the solution point, by assuming them to be linear functions of *x* and *t* inside each individual SE. Since line segments AC and AE belong to (SE) *n*−1/2 *<sup>j</sup>*−1/<sup>2</sup>, *UL* and *FL* can be obtained in terms of the known information stored at (*j* − 1/2, *n* − 1/2), i.e., point A in Fig. 2.8. Performing Taylor expansions in (SE) *n*−1/2 *<sup>j</sup>*−1/<sup>2</sup>yields

$$U\_L = u\_{j-1/2}^{n-1/2} + \frac{\Delta x}{4} (u\_x)\_{j-1/2}^{n-1/2},\tag{2.21}$$

$$F\_L = f\_{j-1/2}^{n-1/2} + \frac{\Delta t}{4} (f\_t)\_{j-1/2}^{n-1/2} \tag{2.22}$$

Similarly, because BC and BF belong to (SE) *n*−1/2 *<sup>j</sup>*+1/<sup>2</sup>, *UR* and *FR* can be obtained by the Taylor expansions about solution point (*j* + 1/2, *n* − 1/2):

$$U\_{\mathcal{R}} = u\_{j+1/2}^{n-1/2} - \frac{\Delta x}{4} (u\_x)\_{j+1/2}^{n-1/2} \tag{2.23}$$

$$F\_{\mathbb{R}} = f\_{j+1/2}^{n-1/2} + \frac{\Delta t}{4} (f\_t)\_{j+1/2}^{n-1/2} \tag{2.24}$$

In Eqs. (2.21)–(2.24), *u* and *ux* at time level *n* − 1/2 are the known values. In addition, since *f* = *f* (*u*) is a prescribed function in the conservation law,

$$f\_{j \pm 1/2}^{n-1/2} = f(u\_{j \pm 1/2}^{n-1/2}) \tag{2.25}$$

By applying the chain rule, the spatial derivative of *f* is described as

16 2 Non-dissipative Core Scheme of CESE Method

$$(f\_{\mathbf{x}})\_{j \pm 1/2}^{n-1/2} = \left[ (\partial f \slash \partial \mathbf{u}) \mu\_{\mathbf{x}} \right]\_{j \pm 1/2}^{n-1/2} \tag{2.26}$$

The temporal derivative of *u* can be obtained using Eq. (2.12),

$$(u\_t)\_{j \pm 1/2}^{n-1/2} = -(f\_x)\_{j \pm 1/2}^{n-1/2} \tag{2.27}$$

By applying the chain rule again, the temporal derivative of *f* is derived as

$$\left( \left( f\_t \right)\_{j \pm 1/2}^{n-1/2} = \left[ \left( \partial f \right/\partial \mu \right) \mu\_t \right]\_{j \pm 1/2}^{n-1/2} \tag{2.28}$$

Accordingly, *UL*, *FL*, *UR*, and *FR* can be explicitly calculated using the above formulas.

The first-order Taylor expansion in (SE) *n <sup>j</sup>*is used again to relate *UL*  \* and *UR*  \* to unknowns *u<sup>n</sup> <sup>j</sup>*and (*ux* ) *n <sup>j</sup>*at time level *n*:

$$U\_L^\* = u\_j'' - \frac{\Delta x}{4} (u\_x)\_j'' \tag{2.29}$$

$$U\_R^\* = u\_j^n + \frac{\Delta x}{4} (u\_x)\_j^n \tag{2.30}$$

Substituting Eqs. (2.29) and (2.30) into (2.19) and (2.20) yields

$$(u\_j'' - \frac{\Delta x}{4}(u\_x)\_j'' = U\_L + (F\_L - F\_C)\frac{\Delta t}{\Delta x},\tag{2.31}$$

$$
\mu\_j'' + \frac{\Delta x}{4} (\mu\_x)\_j'' = U\_R + (F\_C - F\_R) \frac{\Delta t}{\Delta x} \tag{2.32}
$$

Adding Eqs. (2.31)–(2.32), one can derive

$$
\mu\_j'' = \frac{1}{2}(U\_L + U\_R) + \frac{\Delta t}{2\Delta x}(F\_L - F\_R) \tag{2.33}
$$

and subtracting Eqs. (2.31) from (2.32), one has

$$\frac{\Delta x}{4}(u\_x)\_j^n = \frac{1}{2}(U\_R - U\_L) + \frac{\Delta t}{2\Delta x}(2F\_C - F\_L - F\_R) \tag{2.34}$$

Up to now, Eq. (2.33) is presented as an explicit time-marching formula for *u<sup>n</sup> j* , but (*ux* ) *n <sup>j</sup>*still cannot be directly calculated by Eq. (2.34). The reason is because *FC*, which presents the average value of flux *f* through the interface CD, has not been addressed. In the original CESE scheme, *FC* is provided by the Taylor expansion in (SE) *n j* , i.e.,

#### 2.3 Non-dissipative Core Scheme: a Scheme 17

$$F\_C = f\_j'' - \frac{\Delta t}{4} (f\_t)\_j'' \tag{2.35}$$

With a similar procedure to Eqs. (2.25)–(2.28), *f <sup>n</sup> <sup>j</sup>*and ( *ft* ) *n <sup>j</sup>*can be expressed in terms of *u<sup>n</sup> <sup>j</sup>*and (*ux* ) *n j* . Since *u<sup>n</sup> <sup>j</sup>*is determined by Eq. (2.33), a time-marching formula for the only unknown in Eq. (2.34), (*ux* ) *n <sup>j</sup>*, can eventually be derived.

Note that a complete CESE time step consists of two half steps: a half step marching from nodes to centres and another half step marching from centres to nodes. The marching schemes for both half steps are identical, except for the indexes of the solution points. Usually, the initial conditions of *u* and *ux* are specified at the initial time, *t*0 = 0 and the corresponding boundary conditions for *u* and *ux* need to be implemented at the boundaries of *x*.

The above CESE scheme is referred to as the *a* scheme [1] in the literature. Define a solution vector

$$\mathbf{q}\_{j}^{n} = \begin{bmatrix} u\_{j}^{n} \\ \frac{\Delta x}{4} (u\_{x})\_{j}^{n} \end{bmatrix},\tag{2.36}$$

and let function *f* (*u*) be a linear one:

$$f = au, \ a \text{ is a constant} \tag{2.37}$$

The *a* scheme, which is mainly expressed by Eqs. (2.33) and (2.34), can then be rewritten in a matrix form:

$$\mathbf{q}\_{j}^{n} = \mathbf{Q}\_{L}\mathbf{q}\_{j-1/2}^{n-1/2} + \mathbf{Q}\_{R}\mathbf{q}\_{j+1/2}^{n-1/2} \tag{2.38}$$

where the coefficient matrices are

$$\mathbf{Q}\_{L} = \frac{1}{2} \begin{bmatrix} 1+\nu & 1-\nu^{2} \\ -1 & -1+\nu \end{bmatrix}, \mathbf{Q}\_{R} = \frac{1}{2} \begin{bmatrix} 1-\nu & -1+\nu^{2} \\ 1 & -1-\nu \end{bmatrix} \tag{2.39}$$

and ν is a constant, which is defined as

$$\nu \equiv a \frac{\Delta t}{\Delta x} \tag{2.40}$$

Analogously, we apply Eq. (2.38) to the half-step marching from time level *n*−1 to time level *n*−1/2, and obtain

$$\mathbf{q}\_{j-1/2}^{n-1/2} = \mathbf{Q}\_L \mathbf{q}\_{j-1}^{n-1} + \mathbf{Q}\_R \mathbf{q}\_j^{n-1},\tag{2.41}$$

$$\mathbf{q}\_{j+1/2}^{n-1/2} = \mathbf{Q}\_L \mathbf{q}\_j^{n-1} + \mathbf{Q}\_R \mathbf{q}\_{j+1}^{n-1},\tag{2.42}$$

After substituting Eqs. (2.41) and (2.42) into (2.38), a complete time step in the *a* scheme can be formulated as.

$$\mathbf{q}\_{j}^{n} = (\mathbf{Q}\_{L})^{2}\mathbf{q}\_{j-1}^{n-1} + (\mathbf{Q}\_{L}\mathbf{Q}\_{R} + \mathbf{Q}\_{R}\mathbf{Q}\_{L})\mathbf{q}\_{j}^{n-1} + (\mathbf{Q}\_{R})^{2}\mathbf{q}\_{j+1}^{n-1}.\tag{2.43}$$

In fact, the numerical dissipation, dispersion, and stability of different CESE schemes can be conveniently analysed using their matrix forms like Eq. (2.43).

An unique feature of this *a* scheme is its use of Taylor expansion in the inverse time direction, to relate the interface flux *FC* to the marching variables (see Eq. (2.35)). Such a treatment makes the *a* scheme the space–time inversion invariant and renders this scheme non-dissipative [2]. From a historical perspective, the *a* scheme forms the non-dissipative core of subsequent CESE-family schemes.

The *a* scheme described in this chapter is restricted to 1D scalar problems and uniform meshes. First-order Taylor expansions are used in each solution element in the present scheme. Since the *a* scheme is reversible in time, it cannot be used for practical applications, considering that real physical processes are irreversible in time to respect the second law of thermodynamics. Accordingly, substantial improvements have been made to develop this *a* scheme into a viable numerical method. In Chap. 3, we will describe different approaches to introducing necessary numerical dissipation into the CESE scheme, which result in two categories of CESE schemes: the central CESE schemes and the upwind CESE schemes. In Chap. 4, we will describe the extensions of CESE schemes to multi-dimensional Cartesian and unstructured computational meshes. In Chap. 5, high-order Taylor expansions in solution elements will be used to achieve high-order accuracy in both space and time, without enlarging the spatial stencil or adding time-integration stages.

#### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **Chapter 3 CESE Schemes with Numerical Dissipation**

As depicted in Chap. 2, the interface between the two sub-CEs (CD in Fig. 2.7), belongs to the SE of (*j, n*). The flux *FC* needs to be calculated through the Taylor expansion at point (*j, n*) toward the inverse time direction. As a result, the *a* scheme is reversible. This violates the second law of thermodynamics. Thus, the non-dissipative core suffers from the unphysical oscillations for practical applications. This chapter will present the improvements based on the non-dissipative core by introducing necessary numerical dissipation. The resulting dissipative CESE schemes can be categorized into the central CESE schemes and the upwind CESE schemes. It is convenient to begin with the most widely spread CESE scheme, the *a*–α scheme.

#### **3.1** *a–α* **Scheme**

Similar to the *a* scheme described in Chap. 2, the *a*–α scheme is also a central scheme because of the fact that techniques to calculate upwind numerical fluxes are not used. Again, we consider the 1D scalar conservation law of the form

$$\frac{\partial \mu}{\partial t} + \frac{\partial f(u)}{\partial x} = 0 \tag{3.1}$$

To solve Eq. (3.1) with the *a*–α scheme, the discretization and the definitions of CE and SE keep the same as those for the *a* scheme, which can be found in Sect. 2.2. Basically, the solution algorithm follows that of the *a* scheme (see Sect. 2.3), except for Eq. (2.35) for the flux *FC*. It is worth noting that abandoning the calculation of *FC*  does not affect the correctness of Eq. (2.33) for updating *u<sup>n</sup> <sup>j</sup>*. Actually, Eq. (2.33) is purely a result of conservation in the entire conservation element (CE) *n <sup>j</sup>*, irrespective of the flux *FC*. Therefore, the formula for updating the node value of the unknown *u*  remains unchanged, which is written as

<sup>©</sup> The Author(s) 2023 C.-Y. Wen et al., *Space–Time Conservation Element and Solution Element Method*, Engineering Applications of Computational Methods 13, https://doi.org/10.1007/978-981-99-0876-9\_3

22 3 CESE Schemes with Numerical Dissipation

$$
\mu\_j'' = \frac{1}{2}(U\_L + U\_R) + \frac{\Delta t}{2\Delta x}(F\_L - F\_R) \tag{3.2}
$$

The calculation of *UL*, *FL*, *UR*, and *FR* are explicit and has been presented in Sect. 2.3, and then *u<sup>n</sup> <sup>j</sup>*can be evaluated using Eq. (3.2).

The remaining part of the *a*–α scheme is to update the spatial derivative (*ux* ) *n j*  in a different way from the procedure of the *a* scheme. It is in this step where the necessary numerical dissipation is added into the *a*–α scheme. First, two different estimations for (*ux* ) *n <sup>j</sup>*can be obtained:

$$\left(\mu\_x^{-}\right)\_j^n = \frac{u\_j^n - \left[u\_{j-1/2}^{n-1/2} + (\Delta t/2)(u\_t)\_{j-1/2}^{n-1/2}\right]}{\Delta x/2},\tag{3.3}$$

$$\left(\left(u\_{x}^{+}\right)\_{j}^{n} = \frac{\left[u\_{j+1/2}^{n-1/2} + (\Delta t/2)(u\_{t})\_{j+1/2}^{n-1/2}\right] - u\_{j}^{n}}{\Delta x/2} \tag{3.4}$$

In these equations, the temporal derivatives *ut* at time level *n*−1/2 can be readily obtained using Eqs. (2.27) and (2.26). Next, (*ux* ) *n <sup>j</sup>*is taken as a weighted average of (*u*<sup>−</sup> *x* ) *n j* and (*u*<sup>+</sup> *x* ) *n j* :

$$(\boldsymbol{u}\_x)\_j^n = W\left( (\boldsymbol{u}\_x^-)\_j^n, \ (\boldsymbol{u}\_x^+)\_j^n, \ \boldsymbol{\alpha} \right) \tag{3.5}$$

where *W* is a weighted average function with an adjustable parameter α (α ≥ 0, the commonly used values are α = 0, 1 and 2), expressed as

$$W(\mathbf{x}^-, \mathbf{x}^+, \alpha) = \frac{|\mathbf{x}^+|^\alpha \mathbf{x}^- + |\mathbf{x}^-|^\alpha \mathbf{x}^+}{|\mathbf{x}^+|^\alpha + |\mathbf{x}^-|^\alpha},\tag{3.6}$$

The role of this weighted average function is similar to the slope limiter functions in upwind schemes, and it proves to be effective in suppressing the spurious oscillations near a discontinuity.

The above *a*–α scheme has a clear and simple logic. However, the dissipation of the *a*–α scheme increases dramatically as the Courant–Friedrichs–Lewy (CFL) number (the parameter defined by Eq. (2.40), ν = *a*Δ*t*/Δ*x*) approaches zero. In practice, for small values of ν, the discontinuities in the solution can be smeared out. Therefore, the accuracy of the *a*–α scheme is considered to be CFL-number sensitive.

To demonstrate this shortcoming, we consider *f* = *au* (*a* is a constant) in Eq. (3.1), and take α = 0 in Eq. (3.6). In this case, the weighted average function is reduced to a simple arithmetic averaging, and Eq. (3.5) can be written as

$$\left(u\_x\right)\_j^n = \frac{\left[u\_{j+1/2}^{n-1/2} + (\Delta t/2)(u\_t)\_{j+1/2}^{n-1/2}\right] - \left[u\_{j-1/2}^{n-1/2} + (\Delta t/2)(u\_t)\_{j-1/2}^{n-1/2}\right]}{\Delta x},\tag{3.7}$$

which means a central difference approximation at time level *n* is used to update the spatial derivative (*ux* ) *n <sup>j</sup>*. Recall that in Sect. 2.3 we have shown that a CESE scheme applied to the linear scalar convection equation can be written in a matrix form

$$\mathbf{q}\_{j}^{n} = \mathbf{Q}\_{L}\mathbf{q}\_{j-1/2}^{n-1/2} + \mathbf{Q}\_{R}\mathbf{q}\_{j+1/2}^{n-1/2} \tag{3.8}$$

for a half step or

$$\mathbf{q}\_{j}^{n} = (\mathbf{Q}\_{L})^{2}\mathbf{q}\_{j-1}^{n-1} + (\mathbf{Q}\_{L}\mathbf{Q}\_{R} + \mathbf{Q}\_{R}\mathbf{Q}\_{L})\mathbf{q}\_{j}^{n-1} + (\mathbf{Q}\_{R})^{2}\mathbf{q}\_{j+1}^{n-1} \tag{3.9}$$

for a complete step. Here, the solution vector *q* is defined by Eq. (2.36) and the matrices *QL* and *QR* are functions of the CFL number ν = *a*Δ*t*/Δ*x*. By combining Eqs. (3.2) and (3.7), the *a*–α scheme with α = 0 can also be cast into the form of Eq. (3.9) with *QL* and *QR* as follows:

$$\mathbf{Q}\_L = \frac{1}{2} \begin{bmatrix} 1+\nu & 1-\nu^2 \\ -1/2 & \nu \end{bmatrix}, \mathbf{Q}\_R = \frac{1}{2} \begin{bmatrix} 1-\nu & -1+\nu^2 \\ 1/2 & -\nu \end{bmatrix} \tag{3.10}$$

Now, consider a limiting case of Δ*t* = 0 (but Δ*x* > 0), which leads to ν = 0. It can be shown that Eq. (3.9) leads to

$$\mathbf{q}\_{j}^{n} = \frac{1}{4} \begin{bmatrix} 1/2 & 1\\ -1/2 & -1/2 \end{bmatrix} \mathbf{q}\_{j-1}^{n-1} + \frac{1}{4} \begin{bmatrix} 3 & 0\\ 0 & 1 \end{bmatrix} \mathbf{q}\_{j}^{n-1} + \frac{1}{4} \begin{bmatrix} 1/2 & -1\\ 1/2 & -1/2 \end{bmatrix} \mathbf{q}\_{j+1}^{n-1} \tag{3.11}$$

However, a reasonable time-marching scheme should guarantee that **q***<sup>n</sup> <sup>j</sup>*<sup>=</sup> **<sup>q</sup>***n*−<sup>1</sup> *j*  when Δ*t* vanishes, which indicates that the *a*–α scheme suffers from numerical error when the CFL number is very small. Further numerical experiments show that, in practice, when the CFL number ν becomes smaller than 0.1, the excessive numerical dissipation of the *a*–α scheme can be remarkable.

#### **3.2 Courant–Number–Insensitive Scheme**

To overcome the shortcoming of the *a*–α scheme, a Courant–Number–Insensitive (CNI) scheme was constructed [1]. This improvement is based on two requirements: (1) the CNI scheme should reduce to the non-dissipative *a* scheme when the CFL number ν = 0 and (2) the CNI scheme should resemble the *a*–α scheme when |ν| = 1. In this section, the construction of the CNI scheme will be present, and then the relationship between the CNI scheme, the *a* scheme, and the *a*–α scheme will be shown.

First, the formula for updating the node value of the unknown *u* remains the same as Eq. (3.2), which has been used in *a* scheme and the *a*–α scheme. By substituting formulas for *UL*, *FL*, *UR*, and *FR* (Eqs. (2.21)–(2.28)) into Eq. (3.2), the explicit time marching formula for *u<sup>n</sup> <sup>j</sup>*can be written as

$$u\_j^n = \frac{1+\nu}{2} u\_{j-1/2}^{n-1/2} + \frac{1-\nu^2}{8} \Delta x (u\_x)\_{j-1/2}^{n-1/2} + \frac{1-\nu}{2} u\_{j+1/2}^{n-1/2} + \frac{\nu^2 - 1}{8} \Delta x (u\_x)\_{j+1/2}^{n-1/2} \tag{3.12}$$

Then, an algorithm for updating (*ux* ) *n <sup>j</sup>*is devised. As shown in Fig. 3.1, the points (*j* − 1/2, *n*) and (*j* + 1/2, *n*) are denoted by *V*<sup>−</sup> and *V*+, respectively. The point *M*<sup>−</sup> is the midpoint between (*j* − 1/2, *n*) and (*j*, *n*), and *M*+ is the midpoint between (*j*  + 1/2, *n*) and (*j*, *n*). We define two important points *P*<sup>−</sup> and *P*+, which are on the line segments *V* <sup>−</sup>*M*<sup>−</sup> and *V*+*M*+, respectively. The distances between *P*<sup>±</sup> and node (*j*, *n*) are marked in Fig. 3.1. It is clear that *P*<sup>±</sup> coincides with *V*<sup>±</sup> when ν = 1, and they approach *M*<sup>±</sup> as ν → 0. The values of *u*(*x*, *t*) at *P*<sup>±</sup> are estimated by Taylor expansion at points (*j* ± 1/2, *n* − 1/2):

$$
\mu(P\_{\pm}) = \mu\_{j\pm 1/2}^{n-1/2} + \frac{\Delta t}{2} (\mu\_t)\_{j\pm 1/2}^{n-1/2} \mp \frac{(1-|\nu|)\Delta x}{4} (\mu\_x)\_{j\pm 1/2}^{n-1/2} \tag{3.13}
$$

where *un*−1/<sup>2</sup> *<sup>j</sup>*±1/<sup>2</sup>and (*ux* ) *n*−1/2 *<sup>j</sup>*±1/<sup>2</sup>are known, while (*ut* ) *n*−1/2 *<sup>j</sup>*±1/<sup>2</sup>can be obtained using Eqs. (2.27) and (2.26).

Next, *u*(*P*−) and *u*(*P*+) are used to construct two different approximations of (*ux* ) *n j*

$$\left(\hat{\mu}\_x^{-}\right)\_j^n = \frac{\mu\_j^n - \mu(P\_-)}{(1+|\nu|)\Delta x/4} \tag{3.14}$$

and

$$\left(\hat{u}\_x^+\right)\_j^n = \frac{\mu(P\_+) - \mu\_j^n}{(1+|\nu|)\Delta x/4} \tag{3.15}$$

**Fig. 3.1** Definition of points in CNI scheme

#### 3.2 Courant–Number–Insensitive Scheme 25

Finally, the CNI scheme use a weighted average of (*u*ˆ<sup>−</sup> *x* ) *n j* and (*u*ˆ<sup>+</sup> *x* ) *n <sup>j</sup>*as the updated value of (*ux* ) *n <sup>j</sup>*, in a way similar to the *a*–α scheme. However, the weighted average function in the *a*–α scheme, i.e. Eq. (3.6), is replaced with a more sophisticated one as

$$\left(u\_x\right)\_j^n = \frac{\left[1 + f\left(|\nu|\right)(s\_-)\_j^n\right](\hat{u}\_x^+)\_j^n + \left[1 + f\left(|\nu|\right)(s\_+)\_j^n\right](\hat{u}\_x^-)\_j^n}{2 + f\left(|\nu|\right)\left[(s\_-)\_j^n + (s\_+)\_j^n\right]},\tag{3.16}$$

where

$$(s\_{\pm})\_{j}^{n} = \frac{\left| (\hat{u}\_{x}^{\pm})\_{j}^{n} \right|}{\min\left( \left| (\hat{u}\_{x}^{-})\_{j}^{n} \right|, \left| (\hat{u}\_{x}^{+})\_{j}^{n} \right| \right)} - 1 \tag{3.17}$$

and

$$f(|\nu|) = 0.5 \not> |\nu|. \tag{3.18}$$

Hence, the overshoot phenomenon near discontinuities can be suppressed by the artificial dissipation, just like the *a*–α scheme. Moreover, the dissipation brought by Eq. (3.16) can be adjusted dynamically according to the CFL number ν.

When ν = 0, we can show that the CNI scheme will reduce to the non-dissipative *a* scheme. Substituting ν = 0 (which also means Δ*t* = 0) into Eqs. (3.12) and (3.13) yields

$$u\_j^n = \frac{1}{2} u\_{j - 1/2}^{n - 1/2} + \frac{1}{8} \Delta x (u\_x)\_{j - 1/2}^{n - 1/2} + \frac{1}{2} u\_{j + 1/2}^{n - 1/2} - \frac{1}{8} \Delta x (u\_x)\_{j + 1/2}^{n - 1/2} \tag{3.19}$$

and

$$
\mu(P\_{\pm}) = \mu\_{j \pm 1/2}^{n-1/2} \mp \frac{\Delta x}{4} (\mu\_x)\_{j \pm 1/2}^{n-1/2} \tag{3.20}
$$

From Eqs. (3.19) and (3.20), we can get

$$
\mu\_j^n = \frac{\mu(P\_-) + \mu(P\_+)}{2} \tag{3.21}
$$

By using ν = 0 and Eq. (3.21), it is readily shown that

$$\left(\hat{\mu}\_{x}^{+}\right)\_{j}^{n} = \frac{\mu(P\_{+}) - \mu(P\_{-})}{\Delta x/2} = \left(\hat{\mu}\_{x}^{-}\right)\_{j}^{n} \tag{3.22}$$

Thus, the weighted average of (*u*ˆ<sup>−</sup> *x* ) *n j* and (*u*ˆ<sup>+</sup> *x* ) *n <sup>j</sup>*is equal to [*u*(*P*+) −*u*(*P*−)]/(Δ*x*/2), which gives

26 3 CESE Schemes with Numerical Dissipation

$$(u\_x)\_j^n = \frac{1}{\Delta x / 2} \left\{ \left[ u\_{j+1/2}^{n-1/2} - \frac{\Delta x}{4} (u\_x)\_{j+1/2}^{n-1/2} \right] - \left[ u\_{j-1/2}^{n-1/2} + \frac{\Delta x}{4} (u\_x)\_{j-1/2}^{n-1/2} \right] \right\} \tag{3.23}$$

Recall the expression of the *a* scheme (see Eqs. (2.38) and (2.39) in Chap. 2) and let ν = 0, we can reproduce Eqs. (3.19) and (3.23). Therefore, in the limiting case of ν = 0, the CNI scheme and the *a* scheme are identical.

As for the other limiting case |ν| = 1, Eqs. (3.13)–(3.15) lead to

$$\left(\hat{\mu}\_{x}^{-}\right)\_{j}^{n} = \frac{u\_{j}^{n} - \left[u\_{j-1/2}^{n-1/2} + (\Delta t/2)(u\_{t})\_{j-1/2}^{n-1/2}\right]}{\Delta x/2},\tag{3.24}$$

and

$$\left(\hat{\mu}\_{x}^{+}\right)\_{j}^{n} = \frac{\left[\mu\_{j+1/2}^{n-1/2} + (\Delta t/2)(\mu\_{t})\_{j+1/2}^{n-1/2}\right] - \mu\_{j}^{n}}{\Delta x/2} \tag{3.25}$$

There is no difference between (*u*ˆ<sup>±</sup> *x* ) *n j* in the CNI scheme and (*u*<sup>±</sup> *x* ) *n <sup>j</sup>*in the *a*–α scheme. Hence, the algorithms for updating (*ux* ) *n <sup>j</sup>*in the CNI and the *a*–α schemes are basically the same when |ν| = 1, except for the specific forms of the weighted average function.

Both the *a*–α scheme in Sect. 3.1 and the CNI scheme in this section belong to the central CESE schemes. Extensions of these schemes to 2-D and 3-D cases for various systems of conservation equations (e.g. Euler equations for compressible gas dynamics) are straightforward and have been well implemented.

#### **3.3 Upwind CESE Scheme**

The aforementioned central CESE schemes are used to solve nonlinear hyperbolic systems of conservation laws. However, they did not explicitly resort to the knowledge of characteristics or eigenvalues of the systems. As an alternative approach, the characteristic-based upwind CESE scheme proposed by Shen et al. [2], elegantly combines the basic ideas of the CESE method and the upwind numerical flux technique in the Godunov-type FVM method. The upwind CESE scheme is naturally CFL-number-insensitive, which means it does not suffer from the drawback of the *a*–α scheme presented in Sect. 3.1. For some challenging CFD problems such as the simulations of detonations and multiphase flows, the upwind CESE scheme captures discontinuities in flow fields with improved accuracy and robustness, especially for contact discontinuities (e.g., material interfaces).

#### *3.3.1 Construction of Upwind CESE Scheme*

For illustrative purposes, the 1D scalar conservation law (Eq. (3.1)) is considered again. Before introducing the upwind CESE algorithm, some elementary concepts about the space–time discretization need to be clarified. The same computational mesh and solution points as shown in Figs. 2.3 and 2.4 are adopted. The definition of CEs is also retained (see Fig. 2.5). Nevertheless, some modification is made to define the SEs, as sketched in Fig. 3.2. In comparison to the (SE) *n <sup>j</sup>*in Fig. 2.6, the newly designed (SE) *n <sup>j</sup>*no longer contains any part that allows for *t* < *tn*. Apparently, the SEs pave the whole space–time domain without overlap. It is assuemed that *u*(*x*, *t*) and *f* (*x*, *t*) are piecewise linear. For instance, inside (SE) *n <sup>j</sup>*, they can be approximated by the first-order Taylor expansion at point (*j*, *n*). The boundaries of conservation element (CE) *n <sup>j</sup>*belong to three solution elements: DE and DF belong to (SE) *n <sup>j</sup>*, AC and AE belong to (SE) *n*−1/2 *<sup>j</sup>*−1/<sup>2</sup>, and BC and BF belong to (SE) *n*−1/2 *<sup>j</sup>*+1/<sup>2</sup>. The interface between (SE) *n*−1/2 *<sup>j</sup>*−1/<sup>2</sup>and (SE) *n*−1/2 *<sup>j</sup>*+1/<sup>2</sup>, i.e. line segment CD, splits the conservation element (CE) *n j*  into two sub-CEs: (CE−) *n <sup>j</sup>*and (CE+) *n j* .

The construction of the upwind CESE scheme closely follows the framework of the non-dissipative core scheme (i.e., the *a* scheme in Sect. 2.3). Indeed, Eqs. (2.16)–(2.34) still hold for the upwind CESE scheme, and they can be derived by the same procedures as presented in Sect. 2.3. Therefore, just like the *a* scheme, the *a*–α scheme, and the CNI scheme, the formula to update *u<sup>n</sup> <sup>j</sup>*in the upwind CESE scheme remains the same as Eq. (3.2). However, for the purpose of updating (*ux* ) *n j* , the Eq. (2.34) will be utilized, which is the direct result of the space–time integral form of conservation law. For convenience, here we recall this equation:

$$\frac{\Delta x}{4}(u\_x)\_j^n = \frac{1}{2}(U\_R - U\_L) + \frac{\Delta t}{2\Delta x}(2F\_C - F\_L - F\_R) \tag{3.26}$$

Note that this equation is actually discarded in the *a*–α and the CNI schemes, but both the *a* scheme and the upwind CESE scheme make full use of it. Furthermore,

**Fig. 3.2** Definitions of CE and SE for upwind CESE schemes

the key difference between the *a* scheme and the upwind CESE scheme lies in the treatment of *FC* in Eq. (3.26), which denotes the average flux through the interface CD between (CE−) *n <sup>j</sup>*and (CE+) *n <sup>j</sup>*(see Fig. 3.2).

As can be inferred from the new definition of solution elements in Fig. 3.2, the operation used in the *a* scheme to link *FC* with (*ux* ) *n <sup>j</sup>*is forbidden in the upwind CESE scheme. Instead, we intend to obtain *FC* from the known data at time level *n* − 1/2, and then insert *FC* into Eq. (3.26) to get (*ux* ) *n <sup>j</sup>*. Toward this end, a local Riemann problem can be built at the midpoint of CD (denoted by (*j*, *n* − 1/4)), because CD is the interface between (SE) *n*−1/2 *<sup>j</sup>*−1/<sup>2</sup>and (SE) *n*−1/2 *<sup>j</sup>*+1/<sup>2</sup>. To be specific, by Taylor expansion in (SE) *n*−1/2 *<sup>j</sup>*−1/<sup>2</sup>, the left state at point (*j*, *n* − 1/4) is evaluated as

$$\left(\left(u\_L\right)\_j^{n-1/4} = U\_L + \frac{\Delta x}{4} \left(u\_x\right)\_{j-1/2}^{n-1/2} + \frac{\Delta t}{4} \left(u\_t\right)\_{j-1/2}^{n-1/2} \tag{3.27}$$

Meanwhile, Taylor expansion in (SE) *n*−1/2 *<sup>j</sup>*+1/<sup>2</sup>provides the right state at point (*j*, *n* − 1/4):

$$\left(\left(u\_R\right)\_j^{n-1/4} = U\_R - \frac{\Delta x}{4} (u\_x)\_{j+1/2}^{n-1/2} + \frac{\Delta t}{4} (u\_t)\_{j+1/2}^{n-1/2} \tag{3.28}$$

Recall that *UL* and *UR* are given by Eqs. (2.21) and (2.23). According to the definition of solution elements, the values of (*uL* ) *n*−1/4 *<sup>j</sup>* and (*u <sup>R</sup>*) *n*−1/4 *<sup>j</sup>* are not necessarily equal to each other, and in general (*uL* ) *n*−1/4 *<sup>j</sup>* and (*u <sup>R</sup>*) *n*−1/4 *<sup>j</sup>* form the discontinuous initial data on the left and right sides of the "diaphragm" CD.

Next, the average flux *FC* through CD can be evaluated based on the local Riemann problem with initial data (3.27) and (3.28) as

$$F\_C = \hat{F}\left(\left(u\_L\right)\_j^{n-1/4}, \left(u\_R\right)\_j^{n-1/4}\right) \tag{3.29}$$

where *F* ˆ stands for any upwind numerical flux solver, or loosely referred to as Riemann solver. Consequently, the special derivative (*ux* ) *n <sup>j</sup>*can be updated explicitly by Eq. (3.26). It is worth noting that, in the upwind CESE scheme, the time marching formula Eq. (3.2) for *u<sup>n</sup> <sup>j</sup>*has no concern with the upwind procedure above, since flux *FC* never appears in Eq. (3.2).

In case of strong discontinuities, using proper limiters for the derivatives in Eqs. (3.27) and (3.28) bocomes crucial to suppress spurious oscilations. When limiters are used, the reconstruction of (*uL* ) *n*−1/4 *<sup>j</sup>* and (*u <sup>R</sup>*) *n*−1/4 *<sup>j</sup>* is written as

$$\left(u\_L\right)\_j^{n-1/4} = U\_L + \frac{\Delta x}{4}u\_x^L + \frac{\Delta t}{4}u\_t^L,\tag{3.30}$$

$$\left(\left(u\_R\right)\_j^{n-1/4} = U\_R - \frac{\Delta x}{4}u\_x^R + \frac{\Delta t}{4}u\_t^R,\tag{3.31}$$

where *u<sup>L</sup> <sup>x</sup>*, *u <sup>R</sup> x* , *u<sup>L</sup> <sup>t</sup>*, and *u <sup>R</sup> <sup>t</sup>*are the limited derivatives. In this book, the weighted biased averaging procedure limiter (WBAP-L2) [3] is adopted. The limited slopes (spatial derivatives) are

$$\begin{aligned} \mu\_x^L = (\mu\_x)\_{j-1/2}^{n-1/2} W(1, \theta\_1^L, \theta\_2^L), \; \theta\_1^L = \frac{(U\_R - U\_L)/(\Delta x/2)}{(\mu\_x)\_{j-1/2}^{n-1/2}}, \; \theta\_2^L = \frac{(\mu\_x)\_{j+1/2}^{n-1/2}}{(\mu\_x)\_{j-1/2}^{n-1/2}}, \\\\ \mu\_x^R = (\mu\_x)\_{j+1/2}^{n-1/2} W(1, \theta\_1^R, \theta\_2^R), \; \theta\_1^R = \frac{(U\_R - U\_L)/(\Delta x/2)}{(\mu\_x)\_{j+1/2}^{n-1/2}}, \; \theta\_2^R = \frac{(\mu\_x)\_{j-1/2}^{n-1/2}}{(\mu\_x)\_{j+1/2}^{n-1/2}}, \end{aligned} \tag{3.23}$$

where the limiter function is

$$W(1, \theta\_1, \theta\_2) = \begin{cases} \frac{5 + 1/\theta\_1 + 1/\theta\_2}{5 + 1/\theta\_1^2 + 1/\theta\_2^2}, & \text{if } \theta\_1 \text{ and } \theta\_2 > 0\\ 0, & \text{else} \end{cases} \tag{3.34}$$

Once the limited slopes are obtained, the limited temporal derivatives can be calculated by the chain rule and the conservation equation itself, i.e.,

$$
\mu\_t^L = -\mu\_x^L \left. \frac{\partial f}{\partial u} \right|\_{u=U\_L}, \mu\_t^R = -\mu\_x^R \left. \frac{\partial f}{\partial u} \right|\_{u=U\_R} \tag{3.35}
$$

#### *3.3.2 Scheme for Linear Scalar Convection Equation*

Let *f* = *au*, then Eq. (3.1) becomes the linear scalar convection equation. Without loss of generality, *a* is assumed to be a positive constant. For this simple situation, the exact solution of the local Riemann problem with initial data (*uL* ) *n*−1/4 *<sup>j</sup>* and (*u <sup>R</sup>*) *n*−1/4 *j*  can be readily obtained. If the limiter is not used, the interface flux *FC* (Eq. (3.29)) is

$$F\_C = a \left[ U\_L + \frac{\Delta x}{4} (u\_x)\_{j-1/2}^{n-1/2} + \frac{\Delta t}{4} (u\_t)\_{j-1/2}^{n-1/2} \right], (a > 0) \tag{3.36}$$

With this *FC*, the upwind CESE scheme, which is a combination of Eqs. (3.2) and (3.26), can be written in a matrix form as Eqs. (3.8) and (3.9), in which *q* = [*u*, (Δ*x*/4)*ux*] *<sup>T</sup>*and the matrices *QL* and *QR* are functions of the CFL number ν = *a*Δ*t*/Δ*x*:

$$\mathbf{Q}\_L = \frac{1}{2} \begin{bmatrix} 1+\nu & 1-\nu^2 \\ -1+\nu & -1+4\nu-\nu^2 \end{bmatrix}, \mathbf{Q}\_R = \frac{1}{2} \begin{bmatrix} 1-\nu & -1+\nu^2 \\ 1-\nu & -1+\nu^2 \end{bmatrix} \tag{3.37}$$

If the initial data (*uL* ) *n*−1/4 *<sup>j</sup>* and (*u <sup>R</sup>*) *n*−1/4 *<sup>j</sup>* are reconstructed using the limiter, the matrix form shown as Eqs. (3.8) and (3.9) still holds, but the matrices *QL* and *Q<sup>R</sup>* should be

$$\mathbf{Q}\_L = \frac{1}{2} \begin{bmatrix} 1+\nu & 1-\nu^2 \\ -1+\nu & -1+2(1+\phi\_L)\nu + (1-2\phi\_L)\nu^2 \end{bmatrix}, \mathbf{Q}\_R = \frac{1}{2} \begin{bmatrix} 1-\nu & -1+\nu^2 \\ 1-\nu & -1+\nu^2 \end{bmatrix} \tag{3.38}$$

where

$$\phi\_L = W(1, \theta\_1^L, \theta\_2^L), \theta\_1^L = \frac{(U\_R - U\_L)/(\Delta x/2)}{(u\_x)\_{j-1/2}^{n-1/2}}, \theta\_2^L = \frac{(u\_x)\_{j+1/2}^{n-1/2}}{(u\_x)\_{j-1/2}^{n-1/2}}\tag{3.39}$$

The scheme without limiter Eq. (3.37) can be regarded as the special case of Eq. (3.38) with the limiter function φ*L* = 1.

Based on matrices *QL* and *QR* in Eq. (3.38), we can examine the property of the upwind CESE scheme at the limiting case of ν = 0. Recall that the *a*–α scheme becomes very diffusive when <sup>ν</sup> <sup>→</sup> 0 and cannot guarantee that **q***<sup>n</sup> <sup>j</sup>*<sup>→</sup> **<sup>q</sup>***n*−<sup>1</sup> *<sup>j</sup>* as Δ*t* → 0 (but Δ*x* is a finite constant). However, the upwind CESE scheme does not suffer from such a deficiency. A direct evaluation of matrices *QL* and *QR* in Eq. (3.38) as ν → 0 shows

$$\begin{aligned} \mathbf{Q}\_L^2 &\rightarrow \frac{1}{4} \begin{bmatrix} 1 & 1 \\ -1 & -1 \end{bmatrix}^2 = \mathbf{0}, \; \mathbf{Q}\_R^2 \rightarrow \frac{1}{4} \begin{bmatrix} 1 & -1 \\ 1 & -1 \end{bmatrix}^2 = \mathbf{0},\\ \mathbf{Q}\_L \mathbf{Q}\_R + \mathbf{Q}\_R \mathbf{Q}\_L &\rightarrow \frac{1}{4} \begin{bmatrix} 1 & 1 \\ -1 & -1 \end{bmatrix} \begin{bmatrix} 1 & -1 \\ 1 & -1 \end{bmatrix} + \frac{1}{4} \begin{bmatrix} 1 & -1 \\ 1 & -1 \end{bmatrix} \begin{bmatrix} 1 & 1 \\ -1 & -1 \end{bmatrix} = \mathbf{I} \end{aligned} \tag{3.40}$$

Consequently, the time marching scheme for a complete time step (Eq. (3.9)) can always ensure that **q***<sup>n</sup> <sup>j</sup>*<sup>→</sup> **<sup>q</sup>***<sup>n</sup>*−<sup>1</sup> *<sup>j</sup>* as Δ*t* → 0, regardless of the form of the limiter function. Unlike the *a*–α scheme, the upwind CESE scheme is inherently CFLnumber insensitive.

#### *3.3.3 Scheme for Euler Equations*

The 1D unsteady Euler equations for a perfect gas are written as

$$
\frac{
\partial \mathbf{U}
}{
\partial t
} + \frac{
\partial \mathbf{F}
}{
\partial x
} = \mathbf{0},
\tag{3.41}
$$

where *U* = [ρ , ρ*u*, *E*] *<sup>T</sup>*, *F* <sup>=</sup> <sup>|</sup> ρ*u*, ρ*u*<sup>2</sup>+ *p*, (*E* + *p*)*u* |*<sup>T</sup>*are the vectors of the conserved variables and the inviscid flux. Here, ρ, *u*, and *p* represent the density, velocity, and pressure, respectively. The total energy per unit volume is denoted by *E*, which is *E* = *p*/(γ − 1) + ρ*u*2/2 with constant ratio of specific heat γ = 1.4.

The upwind CESE scheme for a scalar equation can be directly applied to solve each component of *U*. However, the upwind procedure to calculate the flux through the interface between (CE−) *n <sup>j</sup>*and (CE+) *n <sup>j</sup>*needs a substantial extension, because of the complexity of the eigen-structure of Eq. (3.41). Still, the flux vector *FC* can be expressed as

$$\mathbf{F}\_{\mathcal{C}} = \hat{\mathbf{F}}\left( (\mathbf{U}\_L)\_j^{n-1/4}, (\mathbf{U}\_R)\_j^{n-1/4} \right), \tag{3.42}$$

where **F** ˆ stands for any appropriate upwind numerical flux solver, e.g., approximate Rimann solvers and other upwind flux functions. Moreover, (**U***L* ) *n*−1/4 *<sup>j</sup>* and (**U***R*) *n*−1/4 *j*  are the data on the two sides of line segment CD in Fig. 3.2, which are reconstructed with an appropriate limiter and serve as the discontinuous initial data of the local Riemann problem. Fortunately, any efficient and robust upwind flux solver [4, 5] exsiting in the FVM literature can be employed to calculate *FC*.

#### *3.3.4 Remarks on Upwind CESE Method*

For the purpose of capturing discontinuities, the upwind CESE method introduces necessary numerical dissipation by the upwind procedure to calculate *FC* (see Eqs. (3.29)–(3.35)). This fact might lead to confusion between the classical upwind FVM [4, 5] and the present upwind CESE method. In the former, the upwind flux technique to tackle local Riemann problems is viewed as the building block of the whole method. In the latter, the upwind flux technique plays a different role.

To shed light on this issue, the interval [*j* − 1/2, *j* + 1/2] is taken as a representative control volume to analyse, which is called "cell" in the FVM, and marked as line segment AB in Fig. 3.2. In this section, three fluxes are related with [*j* − 1/2, *j* + 1/2], denoted by *FL*, *FR*, and *FC*. Apparently, these fluxes can be classified as the flux through the boundary of the control volume (*FL* and *FR*) and the flux through the diaphragm inside the control volume (*FC*). For the upwind CESE scheme, the upwind procedure is performed only to calculate*FC*, while fluxes through the boundaries are not linked to any local Riemann problems. Owing to the timemarching strategy on the staggered mesh, information at points *xj*−1/2 and *xj*+1/2 are ready before evaluating *FL* and *FR*. Thus, the procedure to calculate *FL* and *FR* only involves Taylor expansion inside the SE, the definition of physical flux function *f* (*u*), and the Cauchy-Kowalewski procedure (steps to derive *ut* and *f t* from *ux* as shown in Eqs. (2.26)–(2.28)). On the contrary, the non-staggered FVM employs the upwind flux solver to get all fluxes through boundaries of control volumes, because each *FL*  or *FR* should be treated as the result of a local Riemann problem. Usually, *FC* is not considered in the FVM.

Recall that *FC* is absent in Eq. (3.2). Therefore, the upwind procedure in the CESE scheme never affects the time-marching algorithm for the average value of *u*(*x*, *t*) on [*j* − 1/2, *j* + 1/2], but solely redistributes *u*(*x*, *t*) on [*j* − 1/2, *j* + 1/2] through the calculation of *ux* using Eq. (3.26). However, the upwind FVM relies on two upwind numerical fluxes ˆ*f <sup>j</sup>*−1/2 and ˆ*f <sup>j</sup>*+1/2 to update the average value of *u*(*x*, *t*) on [*j* − 1/2, *j* + 1/2]. This is the reason why we claim that the upwind CESE method utilizes the upwind technique in a different way from the traditional upwind FVM.

#### **3.4 Comparison of Different CESE Schemes**

So far, four CESE schemes, namely the *a* scheme, the *a*–α scheme, the CNI scheme, and the upwind CESE scheme, have been introduced. To make a comparison, their main formulas and properties are listed in the following tables. Table 3.1 indicates that the CNI and upwind CESE schemes are more favourable for the practical simulations, since they can capture discontinuities with reasonable numerical dissipation and they are free of the sensitivity issue which is encountered by the *a*–α scheme. The key formulas in each scheme are shown in Table 3.2. It is notable that all CESE schemes share a common approach to updating *u<sup>n</sup> <sup>j</sup>*. Only the *a* scheme and the upwind CESE scheme update (*ux* ) *n <sup>j</sup>*are based on the conservation law for sub-CEs. In contrast, the *a*–α and the CNI schemes construct (*ux* ) *n <sup>j</sup>*using the finite-difference-like approximation followed by some kind of weighted averaging technique.


**Table 3.1** Properties of four different CESE schemes

**Table 3.2** The formulas to update *u<sup>n</sup> <sup>j</sup>*and (*ux* ) *n <sup>j</sup>*in each CESE scheme



**Table 3.3** The matrices *QL* and *QR* in Eq. (3.8), which is the matrix form of the CESE method applied to the linear scalar convection equation (Eq. (3.1) with *f* = *au*, *a* > 0)

When applied to the linear scalar convection equation (Eq. (3.1) with *f* = *au* and *a* is a positive constant), it is possible to rewrite the CESE scheme in a very compact matrix form. In Table 3.3, the coefficient matrices for different CESE schemes are listed. The matrices for the CNI scheme are not tabulated due to their lengthiness. Such *QL* and *QR*, together with Eqs. (3.8) and (3.9) are useful for revealing the instrinct properties of different CESE schemes, e.g., the stability of each scheme.

#### **3.5 Numerical Examples**

Here, the 1D scalar convection problem and Sod's shock-tube problem with uniform grid size of 0.01 are computed using different CESE schemes. The corresponding C + + source codes implementing the *a*−α CESE scheme are presented in the Appendix. For Eq. (3.1) with *f* = *u*, a square wave propagates in a computational domain [−1, 1] till *t* = 2.0 with the initial condition described as

$$u(x,0) = \begin{cases} 1, & \text{if } -0.5 \le x \le 0.5\\ 0. & \text{else} \end{cases} \tag{3.43}$$

The periodic boundary condition is imposed on both ends. The results are plotted in Fig. 3.3. For CFL = 1, all three schemes provide a virtually exact result. When CFL = 0.8, both the *a*-α and the CNI schemes suffer from strong oscillations on the upwind sides of the discontinuities, while the upwind scheme provides a result with satisfactory accuracy. If CFL is further decreased to an extremely small value, the CNI and upwind scheme can capture the discontinuities very well, but the large dissipation in the *a*−α scheme abnormally smears out the discontinuities. This case proves the CFL insensitivity of the CNI and upwind CESE schemes.

In the Sod's shock-tube problem [6], the gas is initially separated at *x* = 1 with left and right states (ρ , *u*, *p*)*L* = (1.0, 0.0, 1.0) and (ρ , *u*, *p*)*R* = (0.125, 0.0, 0.1). In Fig. 3.4, the computed density profiles at *t* = 0.4 are shown. For CFL = 0.8, all

**Fig. 3.3** Solutions for scalar convection equation: square wave problem at *t* = 2

the schemes can capture the wave structures precisely. For an extremely small CFL number, the CNI and upwind schemes maintain great accuracy, but the *a-*α scheme seriously smeared out the rarefaction, contact, and shock waves.

**Fig. 3.4** Density profile of Sod's shock-tube problem at *t* = 0.4

#### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **Chapter 4 Multi-dimensional CESE Schemes**

The previous chapter has shown that the necessary numerical dissipation can be introduced in 1D CESE schemes through either a central or upwind approach. In this chapter, we present the extensions of 2D CESE schemes based on cartesian meshes, followed by a detailed description of implementation on unstructured meshes. Both the central schemes and upwind schemes will be presented. Then, numerical examples will be provided to demonstrate the capabilities of the present schemes.

#### **4.1 CESE Schemes on Cartesian Meshes**

#### *4.1.1 The Improved a-α CESE Scheme*

The marching scheme requires physical variables and their spatial derivatives at each mesh point. In the pioneer development of the 2D CESE solver [1], the domain is discretized by congruent triangles. Based on a similar technique, Zhang, Yu, and Chang [2] reported a further extension of the 2D CESE scheme on quadrilateral meshes. In the above versions of CESE schemes, the solution points are solely updated at the cell centers with a staggered stencil. As a result, a generalized flux technique was also proposed as a post-marching procedure to handle the "flux decoupling" problem when two mesh points cohost one sub-CE. Alternatively, an improved 2D CESE scheme [3, 4] was proposed with a new definition of SE and CE. In this updated scheme, the entire space–time region is divided into non-overlapping CEs, making it convenient and simpler for calculation and straightforward for extension to 3D scheme. Figure 4.1 shows the space–time geometrical configuration of the improved 2D CESE scheme with a uniform rectangular mesh. The spatial projections of the solution points at (*n* −1/2)th and *n*th time levels are denoted as | and ●, respectively. The solution is updated alternatively between the cell centers ● and cell vertices |. Here Δ*t*/2 = *tn* − *tn*−1/2 and *P*, *P'*, and *P"* subsequently represent the points at three successive timesteps. The conservation element of point *P'*, denoted by CE(*P'*), is defined by the space–time hexahedron *ABCDA'B'C'D'*, i.e., the union of four sub-CEs: CE*LD*(*P'*), CE*LU* (*P'*), CE*RD*(*P'*), and CE*RU* (*P'*). In particular, CE*LD*(*P'*) is defined as *AFPEA'F'P'E'*. Other sub-CEs are defined similarly. Next, consider the solution element of *P'*, SE(*P'*), which consists of three orthogonal planes, *A*'*B*'*C*'*D*', *EGG"E"*, and *FHH"F"*. The physical flux vector is assumed to be smooth within each SE and can be approximated by Taylor expansions about the mesh point associated with the SE. In the following, we present the construction of the 2D CESE solver using similar techniques as described in the 1D scheme.

(a) Mesh points on the *x-y* plane (b) CE(*P*')

**Fig. 4.1** Grid points in the spatial domain, definitions of CE and SE, and corresponding fluxes in a CE for the 2D CESE scheme on rectangular meshes

#### 4.1 CESE Schemes on Cartesian Meshes 39

In this section, we consider the 2D scalar hyperbolic conservation law

$$
\frac{\partial \mu}{\partial t} + \frac{\partial f}{\partial x} + \frac{\partial g}{\partial y} = 0. \tag{4.1}
$$

By using the Gauss' divergence theorem, it is shown that Eq. (4.1) can be written in the form of

$$\oint\_{S(CE(P'))} \mathbf{h} \cdot \mathbf{n} \, \mathrm{d}S = 0,\tag{4.2}$$

where *h* = ( *f*, *g*, *u*) is the space–time flux vector, *S* ( CE( *P*' )) is the boundary of CE( *P*' ) , *n* is the unit outward normal vector on the surface of the control volume. Note that the *u*, *f* , and *g* are approximated by first-order Taylor expansion about point *P'*. For any (*x*, *y*, *t*) ∈ SE( *P*' ) , let

$$
\mu(\mathbf{x}, \mathbf{y}, t) = \mu(\delta \mathbf{x}, \delta \mathbf{y}, \delta t)\_{P'}, \tag{4.3}
$$

$$f(\mathbf{x}, \mathbf{y}, t) = f(\delta \mathbf{x}, \delta \mathbf{y}, \delta t)\_{P'},\tag{4.4}$$

$$\mathbf{g}(\mathbf{x}, \mathbf{y}, t) = \mathbf{g}(\delta \mathbf{x}, \delta \mathbf{y}, \delta t)\_{P'},\tag{4.5}$$

where *X* (δ*x*, δ*y*, δ*t*)*N* denotes the first-order Taylor expansion about point *N* as

$$X(\delta \ge \delta \mathbf{y}, \delta t)\_N = X\_N + (X\_x)\_N \delta \mathbf{x} + \left(X\_\mathbf{y}\right)\_N \delta \mathbf{y} + (X\_\mathbf{t})\_N \delta t,\tag{4.6}$$

and

$$
\delta \mathbf{x} = \mathbf{x} - \mathbf{x}\_N, \ \delta \mathbf{y} = \mathbf{y} - \mathbf{y}\_N, \ \delta t = t - t\_N. \tag{4.7}
$$

In addition, the spatial and temporal derivatives of *f* (*u*) and *g*(*u*) can be derived by the chain rule (Eqs. (2.26) and (2.28)) and *ut* = − *fx* − *gy* .

Meanwhile, one obtains local flux conservation of the four sub-CEs by integrating Eq. (4.2) over their surfaces, given by

$$U\_{LD}^{\prime} \frac{\Delta x \Delta \mathbf{y}}{4} = U\_{LD} \frac{\Delta x \Delta \mathbf{y}}{4} + (F\_{LD} - F\_{CD}) \frac{\Delta \mathbf{y} \Delta t}{4} + (G\_{DL} - G\_{CL}) \frac{\Delta x \Delta t}{4},\tag{4.8}$$

$$U\_{RD}^{\prime} \frac{\Delta x \Delta \mathbf{y}}{4} = U\_{RD} \frac{\Delta x \Delta \mathbf{y}}{4} + (F\_{CD} - F\_{RD}) \frac{\Delta \mathbf{y} \Delta t}{4} + (G\_{DR} - G\_{CR}) \frac{\Delta x \Delta t}{4},\tag{4.9}$$

40 4 Multi-dimensional CESE Schemes

$$U\_{RU}^{\prime} \frac{\Delta x \Delta \mathbf{y}}{4} = U\_{RU} \frac{\Delta x \Delta \mathbf{y}}{4} + (F\_{CU} - F\_{RU}) \frac{\Delta \mathbf{y} \Delta t}{4} + (G\_{CR} - G\_{UR}) \frac{\Delta x \Delta t}{4},\tag{4.10}$$

$$U\_{RU}^{\prime} \quad \Delta x \Delta \mathbf{y} \quad \varepsilon\_{UV} \quad \Delta x \Delta \mathbf{y} \quad \varepsilon\_{UV} \quad \varepsilon\_{UV} \quad \Delta y \Delta t \quad \varepsilon\_{UV} \quad \varepsilon\_{UV} \quad \Delta x \Delta t$$

$$U\_{LU}'\frac{\Delta\lambda\Delta\mathbf{y}}{4} = U\_{LU}\frac{\Delta\lambda\Delta\mathbf{y}}{4} + (F\_{LU} - F\_{CU})\frac{\Delta\mathbf{y}\Delta\mathbf{i}}{4} + (G\_{CL} - G\_{UL})\frac{\Delta\lambda\Delta\mathbf{i}}{4},\tag{4.11}$$

where *ULD*, *U'LD*, *FLD*, *FCD*, *GDL*, and *GCL* are the average fluxes on the surfaces of CE*LD*(*P'*), *AFPE*, *A'F'P'E'*, *AEE'A'*, *FPP'F'*, *AFF'A'*, and *EPP'E'*, respectively. Similar definitions are applied to the other sub-CEs. By definition, the fluxes *ULD*, *FLD*, and *GDL* can be determined via the Taylor expansion in SE(*A*). It is also important to emphasize that the fluxes *FCD* and *GCL* that denote the fluxes across the "inner" surface of two neighboring sub-CEs, however, are not trivially available due to the presence of discontinuities. In fact, it will be shown later that all the fluxes across the interfaces among sub-CEs vanish in the formulation of the time marching scheme of *u*. Apart from this, based on the assumption that the distribution within each SE is linear, the average fluxes on the exterior surfaces of CE*LD*(*P'*) read

$$U\_{LD} = \mu \left(\frac{\Delta x}{4}, \frac{\Delta y}{4}, 0\right)\_A \tag{4.12}$$

$$U\_{LD}' = \mu \left( -\frac{\Delta x}{4}, -\frac{\Delta y}{4}, 0 \right)\_{P'} \tag{4.13}$$

$$G\_{DL} = g\left(\frac{\Delta x}{4}, 0, \frac{\Delta t}{4}\right)\_A \tag{4.14}$$

$$F\_{LD} = \mathbf{g}\left(0, \frac{\Delta \mathbf{y}}{4}, \frac{\Delta t}{4}\right)\_A \tag{4.15}$$

Fluxes are expressed similarly for the other three sub-CEs.

To proceed, by adding Eqs. (4.8)–(4.11), the fluxes through the interfaces between sub-CEs are well balanced. Consequently, the value of the solution point *P'* can be derived as

$$\begin{aligned} u(P') &= \frac{1}{4}(U\_{LD} + U\_{RD} + U\_{RU} + U\_{LU}) + \frac{\Delta t}{4\Delta x}(F\_{LD} - F\_{RD} - F\_{RU} + F\_{LU}) \\ &+ \frac{\Delta t}{4\Delta y}(G\_{LD} + G\_{RD} - G\_{RU} - G\_{LU}) \end{aligned} \tag{4.16}$$

Moreover, by replacing the average fluxes in Eq. (4.16) with the Taylor expansions about their corresponding SEs centers, we can obtain the explicit time marching scheme for *u*(*P'*) as

#### 4.1 CESE Schemes on Cartesian Meshes 41

$$
\mu(P') = \frac{1}{4}\overline{u} + \frac{\Delta t}{4\Delta x}\overline{f} + \frac{\Delta t}{4\Delta y}\overline{g},\tag{4.17}
$$

where

$$\begin{split} \widetilde{u} &= u\left(\frac{\Delta x}{4}, \frac{\Delta y}{4}, 0\right)\_A + u\left(-\frac{\Delta x}{4}, \frac{\Delta y}{4}, 0\right)\_B + u\left(-\frac{\Delta x}{4}, -\frac{\Delta y}{4}, 0\right)\_C + u\left(\frac{\Delta x}{4}, -\frac{\Delta y}{4}, 0\right)\_D, \\ \widetilde{f} &= f\left(0, \frac{\Delta y}{4}, \frac{\Delta t}{4}\right)\_A + f\left(0, -\frac{\Delta y}{4}, \frac{\Delta t}{4}\right)\_D - f\left(0, \frac{\Delta y}{4}, \frac{\Delta t}{4}\right)\_B - f\left(0, -\frac{\Delta y}{4}, \frac{\Delta t}{4}\right)\_C, \\ \widetilde{g} &= g\left(\frac{\Delta x}{4}, 0, \frac{\Delta t}{4}\right)\_A + g\left(-\frac{\Delta x}{4}, 0, \frac{\Delta t}{4}\right)\_B - g\left(-\frac{\Delta x}{4}, 0, \frac{\Delta t}{4}\right)\_C - g\left(\frac{\Delta x}{4}, 0, \frac{\Delta t}{4}\right)\_D. \end{split}$$

With respect to the spatial derivatives, by using the continuous assumptions at points *A', B', C',* and *D'*, the following relations can be formulated

$$\mu\left(-\frac{\Delta x}{2}, -\frac{\Delta \mathbf{y}}{2}, 0\right)\_{p'} = \mu\left(0, 0, \frac{\Delta t}{2}\right)\_A,\tag{4.18}$$

$$\mu\left(\frac{\Delta x}{2}, -\frac{\Delta y}{2}, 0\right)\_{p'} = \mu\left(0, 0, \frac{\Delta t}{2}\right)\_B,\tag{4.19}$$

$$\mu\left(\frac{\Delta x}{2}, \frac{\Delta y}{2}, 0\right)\_{p'} = \mu\left(0, 0, \frac{\Delta t}{2}\right)\_{\mathcal{C}},\tag{4.20}$$

$$\mu\left(-\frac{\Delta x}{2}, \frac{\Delta \mathbf{y}}{2}, 0\right)\_{p'} = \mu\left(0, 0, \frac{\Delta t}{2}\right)\_D. \tag{4.21}$$

After mathematical manipulation of Eqs. (4.19) and (4.20), (4.18) and (4.21), (4.20) and (4.21), (4.18) and (4.19), respectively, the spatial derivatives at solution point *P'* can be explicitly calculated as

$$u\_x^+(P') = \frac{1}{\Delta x} \Big[ u\left(0, 0, \frac{\Delta t}{2}\right)\_B + u\left(0, 0, \frac{\Delta t}{2}\right)\_C - 2u(P') \Big],\tag{4.22}$$

$$
\mu\_x^-(P') = -\frac{1}{\Delta x} \left[ \mu \left( 0, 0, \frac{\Delta t}{2} \right)\_A + \mu \left( 0, 0, \frac{\Delta t}{2} \right)\_D - 2\mu(P') \right], \tag{4.23}
$$

$$
\mu\_{\rm y}^{+}(P') = \frac{1}{\Delta y} \Big[ \mu \Big(0, 0, \frac{\Delta t}{2} \Big)\_{C} + \mu \Big(0, 0, \frac{\Delta t}{2} \Big)\_{D} - 2\mu(P') \Big], \tag{4.24}
$$

$$
\mu\_\mathbf{y}^-(P') = -\frac{1}{\Delta \mathbf{y}} \left[ \mu \left( 0, 0, \frac{\Delta t}{2} \right)\_A + \mu \left( 0, 0, \frac{\Delta t}{2} \right)\_B - 2\mu(P') \right]. \tag{4.25}
$$

One note that, two derivatives were estimated along both spatial directions (denoted by superscripts + and −). Recall that to a weighted averaged function Eq. (3.6) is again adopted to suppress numerical wiggles, such that

$$
\mu\_x(P') = W(\left(\mu\_x^-\right)\_{P'}, \left(\mu\_x^+\right)\_{P'}, \alpha), \tag{4.26}
$$

$$
\mu\_{\mathbf{y}}(P') = W(\left(\mu\_{\mathbf{y}}^{-}\right)\_{P'}, \left(\mu\_{\mathbf{y}}^{+}\right)\_{P'}, \alpha). \tag{4.27}
$$

The Eqs. (4.17) and (4.22)–(4.27) constitute the first half-step marching scheme from cell vertices to cell centers. The second half-step scheme that marches from centers to vertices follows a similar approach. As described above, constructing the improved 2D *a*–α CESE scheme based on rectangular grids is genuinely simple and could be straightforwardly extended to the 3D scheme.

#### *4.1.2 CNI CESE Scheme*

When applying the above 2D *a*–α schemes [1–4] for solving problems with highly non-uniform meshes, the local CFL number may vary significantly from its stability limit of 1 to approaching 0. Namely, the solutions tend to smear out by decreasing the CFL number. The Courant number insensitive (CNI) scheme was proposed to address this issue to alleviate the sensitivity problem by improving the strategy in updating the derivatives [5]. A new point, *Im*, is designated dynamically moving along the line segment connecting the vertices and the centers *M'm* of each sub-CEs (Fig. 4.2). Analogous to the 1D CNI scheme (Sect. 3.2), if we define the local and global Courant numbers as ν and ν0, the interpolated points are designated to approach the centers when ν/ν0 → 0, and approach the cell vertices when ν/ν0 → 1.

For example, the coordinates of *I*' 1 are

$$
\chi(I\_1') = \frac{\nu}{\upsilon\_0} \chi(A') + \left(1 - \frac{\nu}{\upsilon\_0}\right) \chi(M\_1'), \tag{4.28}
$$

$$\mathbf{y}(I\_1') = \frac{\nu}{v\_0}\mathbf{y}(A') + \left(1 - \frac{\nu}{v\_0}\right)\mathbf{y}(M\_1'). \tag{4.29}$$

If we approximate *u*(*I* ' <sup>1</sup>) at *t* = *t*n from the expansions about *A* at *t* = *t*<sup>n</sup>−1/2 and *P'* at *t* = *t*n as

$$\mu(\delta \mathbf{x}\_{+}, \delta \mathbf{y}\_{+}, \mathbf{0})\_{P'} = \mu(I\_{\mathrm{l}}') = \mu\left(\delta \mathbf{x}\_{-}, \delta \mathbf{y}\_{-}, \frac{\Delta t}{2}\right)\_{\mathrm{A}} \tag{4.30}$$

where δ*x*<sup>+</sup> = *x*(*I* ' <sup>1</sup>) − *x*(*P*' ), δ*y*<sup>+</sup> = *y*(*I* ' <sup>1</sup>) − *y*(*P*' ), δ*x*<sup>−</sup> = *x*(*I* ' <sup>1</sup>) − *x*(*A*), and δ*y*<sup>−</sup> = *y*(*I* ' <sup>1</sup>) − *y*(*A*). Similar relations can be derived for each *I* ' *m*. Combining Eq. (4.30)

for *u*(*I* ' <sup>1</sup>) and that for *u*(*I* ' <sup>4</sup>), one can explicitly formulate two independent conditions for calculating the spatial derivatives of *u*(*P'*). Repeating the same procedure for combinations of *u*(*I* ' <sup>2</sup>) and *u*(*I* ' <sup>3</sup>), *u*(*I* ' <sup>1</sup>) and *u*(*I* ' <sup>2</sup>), *u*(*I* ' <sup>3</sup>) and *u*(*I* ' <sup>4</sup>), we arrive at four sets of spatial derivatives. The weighted average function is then used to calculate the optimal derivatives for capturing the discontinuities.

#### *4.1.3 Upwind CESE Scheme*

In the 2D upwind CESE scheme that Shen et al. [6] proposed, the CE retains the same form as the central CESE scheme. The time marching scheme of *u* is the same as in the *a-*α scheme. The upwind procedure only affects the calculation of the spatial derivatives, and global space–time conservation is ensured. In analogue to the definition in the 1D upwind scheme, for example, SE(*P'*) is defined as cuboid *A'B'C'D'A"B"C"D"* (Fig. 4.3), where the physical flux vector is approximated by the Taylor expansion about point *P'*. The average values at the top surfaces of the four sub-CEs of CE(*P'*) read

$$U\_{LD}' = \mu \left( -\frac{\Delta x}{4}, -\frac{\Delta y}{4}, 0 \right)\_{P'} , \tag{4.31}$$

$$U\_{RD}' = \mu \left(\frac{\Delta x}{4}, -\frac{\Delta y}{4}, 0\right)\_{P'} \tag{4.32}$$

$$U\_{RU}' = \mu \left(\frac{\Delta x}{4}, \frac{\Delta y}{4}, 0\right)\_{P'} , \tag{4.33}$$

**Fig. 4.3** Definition of SE of the upwind CESE scheme

$$U\_{LU}' = \mu \left( -\frac{\Delta x}{4}, \frac{\Delta y}{4}, 0 \right)\_{P'} . \tag{4.34}$$

Referring to Fig. 4.1d, by substituting Eqs. (4.31) in (4.8), (4.32) in (4.9) and subtracting the two, one obtains

$$\mu\_X^D(P') = \frac{2}{\Delta x} \left[ (U\_{RD} - U\_{LD}) + \frac{\Delta t}{\Delta x} \left( 2F\_{CD} - F\_{LD} - F\_{RD} \right) + \frac{\Delta t}{\Delta y} \left( G\_{DR} - G\_{CR} - G\_{DL} + G\_{CL} \right) \right]. \tag{4.35}$$

Equation (4.35) represents the spatial derivative in the *x* direction based on the information of the lower pair of CEs. Note that in contrast to central schemes, the upwind scheme requires additional estimations of inner fluxes (e.g., *FCD*, *GCR*, and *GCL* in Eq. (4.35)) to compute spatial derivatives. Analogously, with the aid of Eqs. (4.31) and (4.34) and by combining with Eqs. (4.10) and (4.11), Eqs. (4.8) and (4.11), Eqs. (4.9) and (4.10), respectively. The following equations for spatial derivatives can be formulated

$$\ln u\_X^U(P') = \frac{2}{\Delta x} \left[ \left( U\_{RU} - U\_{LU} \right) + \frac{\Delta t}{\Delta x} \left( 2F\_{CU} - F\_{LU} - F\_{RU} \right) + \frac{\Delta t}{\Delta y} \left( G\_{CR} - G\_{UR} - G\_{CL} + G\_{UL} \right) \right]. \tag{4.36}$$

$$u\_{\mathcal{Y}}^L(P') = \frac{2}{\Delta y} \left[ (U\_{LU} - U\_{LD}) + \frac{\Delta t}{\Delta x} (F\_{LU} - F\_{CU} - F\_{LD} + F\_{CD}) + \frac{\Delta t}{\Delta y} (2G\_{CL} - G\_{DL} - G\_{UL}) \right]. \tag{4.37}$$

$$u\_{\bar{y}}^R(P') = \frac{2}{\Delta y} \left[ \left( U\_{RU} - U\_{RD} \right) + \frac{\Delta t}{\Delta x} \left( F\_{CU} - F\_{RU} - F\_{CD} + F\_{RD} \right) + \frac{\Delta t}{\Delta y} \left( 2G\_{CR} - G\_{DR} - G\_{UR} \right) \right]. \tag{4.38}$$

The final spatial derivatives are calculated by using WBAP limiter to remove small noises near the strong discontinuities, e.g.,

$$u\_{\lambda}(\boldsymbol{P}') = u\_{\lambda}^{\mathbb{C}}(\boldsymbol{P}')WBAP^{L2}(\mathbb{I},\ \theta\_1,\ \theta\_2) = \frac{1}{2} \Big[ u\_{\lambda}^{D}(\boldsymbol{P}') + u\_{\lambda}^{U}(\boldsymbol{P}') \Big] WBAP^{L2}(\mathbb{I},\ \theta\_1,\ \theta\_2),\tag{4.39}$$

where

$$\begin{aligned} \theta\_1 &= \mu\_x^D(P') / \mu\_x^C(P'), \\\\ \theta\_2 &= \mu\_x^U(P') / \mu\_x^C(P'). \end{aligned}$$

In the expressions for spatial derivatives in Eqs. (4.35)–(4.38), the fluxes on CE(*P'*) surfaces can be easily approximated by the Taylor expansions of the values from *t*  = *tn*−1/2. To complete Eqs. (4.35)–(4.38), the fluxes through the inner boundaries of CE(*P'*), *FCD*, *FCU* , *GCL*, and *GCR*, need to be solved. However, the points at the two sides of the interface between two sub-CEs belong to different SEs. In that case, the fluxes across these interfaces may be discontinuous. For this reason, they can be calculated by an upwind procedure. For example,

$$F\_{CD} = \hat{F}(\mu\_{c0}^-, \mu\_{c0}^+),\tag{4.40}$$

where *F* Λ refers to any efficient Riemann solver, *u*<sup>−</sup> *CD* and *u*<sup>+</sup> *CD* denote the values at the centroid of the interface *FPP'F'*, and they are respectively approximated by

$$
\mu\_{c0}^- = U\_{LD} + \left(\mu\_x^\*\right)\_{LD} \frac{\Delta x}{4} + \left(\mu\_t^\*\right)\_{LD} \frac{\Delta t}{4} \tag{4.41}
$$

$$
\mu\_{c0}^{+} = U\_{RD} - \left(\mu\_{x}^{\*}\right)\_{RD} \frac{\Delta x}{4} + \left(\mu\_{t}^{\*}\right)\_{RD} \frac{\Delta t}{4} \tag{4.42}
$$

As an analogue to the 1D upwind CESE scheme, reconstructed slopes are employed in the calculation of the upwind fluxes to eliminate the spurious oscillations:

$$\left(\mu\_{\boldsymbol{x}}^{\*}\right)\_{LD} = \mu\_{\boldsymbol{x}}(\boldsymbol{A}) \boldsymbol{W} \boldsymbol{B} \boldsymbol{A} \boldsymbol{P}^{L2} \left(\boldsymbol{1}, \ \theta\_{LD}^{1}, \ \theta\_{LD}^{2}\right) \tag{4.43}$$

$$\left(\mu\_x^\*\right)\_{RD} = \mu\_x(B) W B A P^{L2} \left(1, \,\, \theta\_{RD}^1, \,\, \theta\_{RD}^2\right) \tag{4.44}$$

where the parameters in WBAP limiter are

$$\begin{split} \theta\_{LD}^{\mathrm{l}} &= \frac{(U\_{RD} - U\_{LD})/(\Delta x/2)}{\mu\_{\mathrm{x}}(A)}, \quad \theta\_{LD}^{2} = \frac{\mu\_{\mathrm{x}}(B)}{\mu\_{\mathrm{x}}(A)},\\ \theta\_{RD}^{\mathrm{l}} &= \frac{(U\_{RD} - U\_{LD})/(\Delta x/2)}{\mu\_{\mathrm{x}}(B)}, \quad \theta\_{RD}^{2} = \frac{\mu\_{\mathrm{x}}(A)}{\mu\_{\mathrm{x}}(B)}. \end{split} \tag{4.45}$$

The temporal derivatives are hereafter calculated by the chain rule and Eq. (4.1),

$$\left(\left(u\_t^\*\right)\_{LD} = -f\_x\left(U\_{LD}, \left(u\_x^\*\right)\_{LD}\right) - g\_y\left(U\_{LD}, u\_y(A)\right)\right.\tag{4.46}$$

$$\left(\mu\_{\,\,t}^{\*}\right)\_{RD} = -f\_{\,\,x}\left(U\_{RD}, \left(\mu\_{\,\,x}^{\*}\right)\_{RD}\right) - g\_{\,\,y}\left(U\_{RD}, \mu\_{\,\,y}(B)\right) \tag{4.47}$$

#### **4.2 CESE Schemes on Unstructured Meshes**

In order to extend the capability of the CESE scheme to solve problems with complex spatial geometries, CESE schemes on unstructured meshes have been proposed and developed. Notably, both central and upwind versions have been developed for hybrid meshes of triangular and quadrilateral elements. Hereafter, the second order CESE scheme on hybrid meshes [7] is presented here.

#### *4.2.1 a-α CESE Scheme*

Consider a 2D mesh consisting of quadrilateral and triangular elements as sketched in Fig. 4.4. The *Vi* denotes the vertices of the cells. Initially, the value of *u*, and its spatial derivatives, *ux* and *uy*, are provided at these vertices. For each cell, *Cm*  represents its centroid, and *Om* is the mid-point of the cell edge. Meanwhile, *Vi* , *V*' *i* , and *V*'' *<sup>i</sup>*denote points at time step *tn*−1/2, *tn*, *tn*+1/2, respectively. Other space–time points are defined similarly.

The staggered marching strategy is again adopted here. The scheme marches from vertices *V*' *<sup>i</sup>*to cell centroids *C*'' *<sup>m</sup>*from *t* = *tn* to *t* = *tn*+1/2. The corresponding CE for*C*'' *m*  is either a triangular prism or a hexahedron (Fig. 4.4c and d). SE( *C*'' *m* ) ) is defined as four (e.g., SE( *C*'' 1 ) ) or five (e.g., SE(*C*'' 2)) planes intersecting at the centers. Similarly, in the other half marching step, from *t* = *tn*−1/2 to *t* = *tn*, *Cm* and *Om* form a polygonal element of which the centroid is denoted as *Gi*, and the scheme marches from *Cm*  to the *Gi*. Note that *Gi* does not necessarily coincide with *Vi*. Thus, interpolation is required at each marching step. SE(*V*' *<sup>i</sup>*) is defined as *M* vertical planes and one horizontal plane intersecting at *V*' *<sup>i</sup>*, where *M* denotes the number of *Cm* linked to CE(*Vi*). Each CE is composed of a group of hexahedral sub-CEs. As shown in Fig. 4.5, *Om*−<sup>1</sup>*Cm Om Vi O*' *<sup>m</sup>*−<sup>1</sup>*C*' *<sup>m</sup> O*' *<sup>m</sup> V*' *i* is the *m*th sub-CE belonging to CE(*V*' *<sup>i</sup>*), *G*˜ *<sup>m</sup>* and *G*˜ ' *<sup>m</sup>*are defined as the centroids of the bottom surface *Om*−<sup>1</sup>*Cm Om Vi* and top surface *O*' *<sup>m</sup>*−<sup>1</sup>*C*' *<sup>m</sup> O*' *<sup>m</sup> V*' *i* , respectively.

**Fig. 4.4** Definition of conservation elements in CESE schemes based on a hybrid mesh

It should be noted that the marching scheme from *t* = *tn* to *t* = *tn*+1/2 is a particular case of that from *t* = *tn*−1/2 to *t* = *tn*. From *t* = *tn* to *t* = *tn*+1/2, the projection of a CE on the spatial domain is either a quadrilateral or a triangle, while the projection in the other step is a polygon with an arbitrary shape. Without loss of generality, we present the marching scheme from *t* = *tn*−1/2 to *t* = *tn*. By applying the integral conservation law Eq. (4.1) to the *m*th sub-CE, the flux-balancing relation is obtained as

$$\iint\limits\_{O\_{n-1}^{\prime}\subset C\_{n}\mathcal{O}\_{n}^{\prime}V\_{i}^{\prime}} u \, d\sigma = \iint\limits\_{O\_{n-1}\subset C\_{n}\mathcal{O}\_{n}\mathcal{V}\_{i}} u \, d\sigma - \iint\limits\_{\begin{subarray}{c}\mathcal{O}\_{n-1}\subset C\_{n}\mathcal{O}\_{n}^{\prime}\mathcal{O}\_{n-1}^{\prime} + \mathcal{C}\_{n}\mathcal{O}\_{n}\mathcal{O}\_{n}^{\prime}\mathcal{C}\_{n}^{\prime}\\ + \mathcal{O}\_{n-1}V\_{i}V\_{i}^{\prime}O\_{n-1}^{\prime} + \mathcal{O}\_{n}V\_{i}V\_{i}^{\prime}O\_{n}^{\prime}\end{subarray}} \mathbf{V} \cdot \mathbf{n} \, d\sigma,\quad(4.48)$$

**Fig. 4.5** The *m*th sub-CE of CE(*Vi*)

where **V** = ( *f*, *g*) is the vector of the fluxes, and **n** is the unit outward normal vector of the surface. In second-order schemes, the integration over a surface is equivalent to the average value at the centroid multiplied by the surface area. Therefore, the above equation can be rewritten as

$$U\_m' = U\_m - \frac{\Delta t}{2S\_m} \sum\_{l=1}^2 \left( \Delta \mathbf{y}\_m^{(l)} F\_m^{(l)} - \Delta x\_m^{(l)} G\_m^{(l)} \right) + \frac{\Delta t}{2S\_m} \left( L\_{m-1} \widehat{F}\_{m-1} - L\_m \widehat{F}\_m \right), \tag{4.49}$$

where *Sm* denotes the area of the top surface, and Δ*t*/2 = *tn* − *tn*−1/2. *U*' *<sup>m</sup>*, *Um*, (*F*(1) *m* , *G*(1) *m* ), and (*F*(2) *m* , *G*(2) *m* ) represent the average values over *O*' *<sup>m</sup>*−<sup>1</sup>*C*' *<sup>m</sup> O*' *<sup>m</sup> V*' *i* , *Om*−<sup>1</sup>*Cm Om Vi* , *Om*−<sup>1</sup>*CmC*' *<sup>m</sup> O*' *<sup>m</sup>*−1, and *Cm Om <sup>O</sup>*' *mC*' *<sup>m</sup>*, respectively. They can be approximated by first-order Taylor expansion in the corresponding SEs. *F* Λ *<sup>m</sup>*−1 and *F* Λ *<sup>m</sup>*are the fluxes in the normal direction across the interfaces *Om*−<sup>1</sup>*ViV*' *i O*' *<sup>m</sup>*−1 and *Om ViV*' *i O*' *<sup>m</sup>*with the neighbouring sub-CEs. *Lm*-1 and *Lm* are the lengths of the line segments *Om*−<sup>1</sup>*Vi* and *Om Vi* , respectively. The summation of the flux-balancing equations of every sub-CEs is expressed as

$$\iint\_{C\_1'O\_2'C\_2''...C\_M'} u \, d\sigma = \sum\_{m=1}^M \iint\_{O\_{m-1}'C\_m'O\_m'V\_i'} u \, d\sigma \tag{4.50}$$

where *O*0 ≡ *OM* . By substituting Eqs. (4.48) and (4.49) into (4.50), the interface fluxes between sub-CEs are cancelled, the conservation law on the entire CE(*V*' *i* ) yields

$$\mu(G\_i') = \left[\sum\_{m=1}^{M} U\_m S\_m - \frac{\Delta t}{2} \sum\_{m=1}^{M} \sum\_{l=1}^{2} \left(\Delta \mathbf{y}\_m^{(l)} F\_m^{(l)} - \Delta \mathbf{x}\_m^{(l)} G\_m^{(l)}\right)\right] \bigg/ \sum\_{m=1}^{M} S\_m \quad (4.51)$$

Equation (4.51) is the explicit marching scheme for *u*(*G*' *i* ). In addition, the treatment of spatial derivatives is similar to that in cartesian meshes. By using the continuous assumption at *C*' *<sup>m</sup>*, the following relation can be established

$$
\mu(C\_m) + \mu\_l(C\_m)\frac{\Delta t}{2} = \mu(C\_m') = \mu(G\_i') + \mu\_x(G\_i')\delta x\_m + \mu\_y(G\_i')\delta y\_m \tag{4.52}
$$

for each *m*, where δ*xm* = *x*(*Cm*) − *x*(*Gi*), and δ*ym* = *y*(*Cm*) − *y*(*Gi*). Apply the same relation at *C*' *<sup>m</sup>*+1, together with Eq. (4.52), the two spatial derivatives are ready to be solved using Cramer's rule,

$$\mu\_x(G\_i') = \frac{\begin{vmatrix} \delta \boldsymbol{u}\_m & \delta \mathbf{y}\_m \\ \delta \boldsymbol{u}\_{m+1} \; \delta \mathbf{y}\_{m+1} \\ \hline \\ \delta \boldsymbol{x}\_{m+1} & \delta \mathbf{y}\_{m+1} \end{vmatrix}}{\begin{vmatrix} \delta \boldsymbol{x}\_m & \delta \mathbf{y}\_m \\ \delta \boldsymbol{x}\_{m+1} & \delta \mathbf{y}\_{m+1} \end{vmatrix}}, \qquad \mu\_y(G\_i') = \frac{\begin{vmatrix} \delta \boldsymbol{x}\_m & \delta \boldsymbol{u}\_m \\ \delta \boldsymbol{x}\_{m+1} & \delta \mathbf{y}\_{m+1} \\ \hline \\ \delta \boldsymbol{x}\_{m+1} & \delta \mathbf{y}\_{m+1} \end{vmatrix}}{\begin{vmatrix} \delta \boldsymbol{x}\_m & \delta \mathbf{y}\_m \\ \delta \boldsymbol{x}\_{m+1} & \delta \mathbf{y}\_{m+1} \end{vmatrix}},\tag{4.53}$$

where

$$
\delta u\_m = \mu(C\_m) + \mu\_t(C\_m)\frac{\Delta t}{2} - \mu(G\_i'). \tag{4.54}
$$

Totally, *M* pairs of derivatives can be calculated in this approach. Shen et al. [7] suggested using the weighted average function to obtain the final derivatives as

$$\mu\_x(G\_i') = \frac{\sum\_{m=1}^M W^{(m)} u\_x^{(m)}(G')}{\sum\_{m=1}^M W^{(m)}}, \quad \mu\_y(G\_i') = \frac{\sum\_{m=1}^M W^{(m)} u\_y^{(m)}(G')}{\sum\_{m=1}^M W^{(m)}},\tag{4.55}$$

where

$$W^{(m)} = \left(\prod\_{\substack{i=1,\ i \neq m}}^M \sqrt{\left[u\_x^{(i)}(G')\_i\right]^2 + \left[u\_y^{(i)}(G')\_i\right]^2}\right)^\alpha \tag{4.56}$$

Here, the adjustable variable α (α = 0 ∼ 2) controls the dissipation near the discontinuities. Smaller α results in lower dissipation, and the function will reduce to an arithmetic average when α = 0. Finally, the value at *V*' *i* is interpolated from *G*' *<sup>i</sup>*as

$$\mu(V\_i') = \mu(G\_i') + \mu\_x(G\_i') \left[ \mathbf{x}(V\_i') - \mathbf{x}(G\_i') \right] + \mu\_y(G\_i') \left[ \mathbf{y}(V\_i') - \mathbf{y}(G\_i') \right]. \tag{4.57}$$

The *a*–α CESE scheme on the hybrid mesh formulated by Eqs. (4.51) and (4.57), while robust, is nevertheless sensitive to the Courant number. The following section addresses the further improvements to make the scheme Courant number insensitive, including the CNI scheme and upwind scheme. These two updated schemes, being closely analogous to those presented in Sects. 4.1.2 and 4.1.3, possess the identical time-marching scheme for *u* as provided in the *a*-α CESE scheme on a hybrid mesh. Only the marching scheme for spatial derivatives is modified accordingly.

#### *4.2.2 CNI CESE Scheme*

As a straightforward extension of the 1D CNI scheme, Shen and Parsani [8] extended the CNI scheme to a hybrid mesh. Recall that the central idea of the CNI scheme is that when the local Courant number (ν) approaches the global Courant number (ν0), the scheme becomes the *a*-α scheme. When ν approaches 0, the scheme reduces to the non-dissipative counterpart. As depicted in Fig. 4.6, a new point *I* ' *<sup>m</sup>*is defined such that

$$\chi(I'\_m) = \frac{\nu}{\upsilon\_0} \chi(C'\_m) + \left(1 - \frac{\nu}{\upsilon\_0}\right) \chi(\tilde{G}'\_m),\tag{4.58}$$

$$\mathbf{y}(I\_m') = \frac{\nu}{\upsilon\_0} \mathbf{y}(C\_m') + \left(1 - \frac{\nu}{\upsilon\_0}\right) \mathbf{y}(\tilde{G}'\_m). \tag{4.59}$$

The value of *I* ' *m* can be obtained using the expansion from *G*' *<sup>i</sup>*as

$$\mu(I\_m') = \mu\left(\mathbf{x}(I\_m') - \mathbf{x}(G\_i'), \mathbf{y}(I\_m') - \mathbf{y}(G\_i'), \mathbf{0}\right)\_{G\_i'},\tag{4.60}$$

**Fig. 4.6** Definition of points in instructed mesh for CNI scheme

where *u*(*I* ' *<sup>m</sup>*) can be calculated as

$$
\mu(I\_m') = \mu\left(\mathbf{x}(I\_m') - \mathbf{x}(C\_m), \mathbf{y}(I\_m') - \mathbf{y}(C\_m), \frac{\Delta t}{2}\right)\_{C\_m}.\tag{4.61}
$$

With the aid of Eqs. (4.60) and (4.61) to each *I* ' *<sup>m</sup>*. Every two neighboring *I* ' *m*  explicitly provide one set of *ux* (*G*' *<sup>i</sup>*) and *uy* (*G*' *<sup>i</sup>*). Then, the optimal derivatives can be weightedly averaged from these *M* sets of derivatives.

#### *4.2.3 Upwind CESE Scheme*

In the upwind CESE scheme, the definition of SE is modified and SE(*Vi*) is the polyhedron *C*' 1*O*' 1*C*' 2...*C*' *<sup>M</sup> O*' *<sup>M</sup> C*'' 1 *P*''*O*'' 1*C*'' 2...*C*'' *<sup>M</sup> O*'' *<sup>M</sup>*(Fig. 4.4b). The values of *u*, *f* , and *g* are continuous within each SE. However, the interface flux *F* Λ *<sup>m</sup>*between two neighboring sub-CEs may be discontinuous. The flux *F* Λ *<sup>m</sup>*can be evaluated in an upwind manner in the normal direction as

$$
\hat{F}\_m = \hat{F}(\mu\_L, \mu\_R) \tag{4.62}
$$

where *F* ˆ denotes the Riemann solver, and the values at the center of the interface, *uL* and *uR*, are reconstructed from SE(*Cm*) and SE(*Cm*+1) with a similar procedure as the 2D upwind CESE scheme on a uniform mesh. The WBAP limiter (Eq. (3.42)) is used to reconstruct the derivatives. Once the fluxes across inner interfaces are determined, the averaged value on the top surface of *m*th sub-CE *U*' *<sup>m</sup>*in Eq. (4.49) can be explicitly computed. We can therefore relate *U*' *m* to *u*(*G*' *<sup>i</sup>*) using the Taylor expansion as

$$\left[U\_m' = \mu(G\_i') + \mu\_x(G\_i')\right] \mathbf{x}(G\_m') - \mathbf{x}(G\_i') \left[ + \mu\_y(G\_i') \left[ \mathbf{y}(G\_m') - \mathbf{y}(G\_i') \right] \right] \tag{4.63}$$

Two linear equations can be derived for each set of *m* and *m* + 1. Accordingly, a similar procedure for averaging the derivatives in the central scheme can be directly applied here.

#### **4.3 Numerical Examples**

Two canonical problems are selected to demonstrate the performance of the *a*–α, CNI, and upwind CESE schemes on either structured or unstructured meshes. For Euler equations, the above schemes can be directly applied to solve each component of the conserved variables. The interface fluxes in the upwind CESE scheme can be solved as a Riemann problem, where Riemann solvers including Harten, Lax and van Leer (HLL) Riemann solver, the contact discontinuity restoring HLLC Riemann solver, and the Roe Riemann solver, can be applied to solve the local Riemann problem. Since the upwind direction may not be perpendicular to the grid-aligned direction, here we use the rotated HLLC Riemann solver. The first example is the Mach 3 wind tunnel problem with step. The computational domain is [0, 3] × [0, 1] with a step of height 0.2 placed at 0.6 from the entrance. The inflow boundary condition is imposed at the left boundary, and supersonic outflow is applied at the right boundary. The upper and lower boundaries are treated as reflective boundaries. Initially, the primitive variables (ρ , *u*, v, *p*) are set as (1.4, 3, 0, 1). Figure 4.7 depicted the simulations using different CESE schemes with 1800 × 600 structured cartesian meshes. Both the CNI scheme and the upwind scheme can capture the shear layer instability that is missed in the result from *a*–α scheme.

The second example is the double Mach reflection problem. A Mach 10 shock propagates to the right and reflects over a solid wedge with an incline angle of 30 ◦. A supersonic inflow boundary condition was implemented on the left boundary. Reflective and non-reflective boundary conditions were applied on the lower and upper boundary, respectively. The problem was solved using unstructured quadrilateral meshes with a mesh size of 1/400. Figure 4.8 compares the results by applying different schemes at computation time *t* = 0.2. All the three schemes can accurately describe the structures while the upwind scheme surpasses with finer resolved feature near the wall-adjacent jet.

**Fig. 4.7** Mach 3 wind tunnel with a step simulated by different CESE schemes with 1800 × 600 mesh grids: density contours at *t* = 4.0

#### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **Chapter 5 High-Order CESE Schemes**

This chapter is dedicated to the description of high-order CESE schemes. In the second-order CESE schemes, the first-order Taylor expansion was employed to approximate the unknowns and fluxes within the solution elements. The accuracy of the scheme mainly depends on the approximations on the surfaces of the conservation elements. Analogously, the high-order CESE schemes with *M*th-order in space and time are generally derived from (*M*−1)th-order Taylor expansions in the solution elements. The high-order CESE schemes use a highly compact stencil. Furthermore, spatial and temporal high-order accuracy can be achieved simultaneously. We shall start with constructing a 1D high-order scheme and then extend it to multi-dimensional schemes.

#### **5.1 Construction of High-Order CESE Schemes**

Several attempts have been made to obtain higher-order accuracy for CESE schemes. One of the primary advantages of the high-order CESE schemes is the usage of the most compact stencil. High-order accuracy is achieved by the approximation with high-order Taylor expansions. For example, *a*(3) scheme [1] has 4th order of accuracy by applying second-order Taylor expansion, and it is stable for *v* < 0.5. Furthermore, *a*(4) scheme [2] with 4th to 5th order of accuracy was developed by defining more CEs at each grid point. The CFL number needs to be constrained below 1/3. The constructions of high-order schemes from the above two approaches were limited to only 1D scenarios. Alternatively, Chang [3] proposed a high-order scheme that can be extended to arbitrary order with CFL limited below 1. In this scheme, the even derivatives are advanced in the same manner as the conserved variables, whereas the odd derivatives are advanced through finite difference. This method has been extended to solve Euler equations [4]. In this section, we follow the approaches of Liu and Wang [5] in deriving the high-order 1D scheme, Wang et al. [6] for 2D high-order schemes on uniform meshes, and Shen et al. [7] for the 2D high-order scheme for hybrid meshes. In these schemes, only the conserved variables are computed through the integration on the CE surfaces, while all the physical derivatives are computed through the finite difference of lower-order derivatives. In the following text, we will present the constructions of the third-order CESE schemes. The construction of arbitrary order schemes could be achieved by implementing the same strategy, and details can be accessed in Shen et al. [7] and Yang et al. [8].

#### *5.1.1 Construction of a Third-Order 1D CESE Scheme*

Recall the 1D scalar problem

$$\frac{\partial \mu}{\partial t} + \frac{\partial f(u)}{\partial x} = 0 \tag{5.1}$$

with *f* = *au*. To formulate a third-order scheme, *u*(*x*, *t*) and *f* (*x*, *t*) in any (*x*, *t*) ∈ SE(*j*, *n*) are approximated by the second-order Taylor expansion as

$$\begin{aligned} \mu(\mathbf{x},t) &= u\_j^n + (u\_x)\_j^n(\mathbf{x}-\mathbf{x}\_j) + (u\_t)\_j^n(t-t\_n) + (u\_{xt})\_j^n(\mathbf{x}-\mathbf{x}\_j)(t-t\_n) \\ &+ \frac{1}{2}(u\_{xx})\_j^n(\mathbf{x}-\mathbf{x}\_j)^2 + \frac{1}{2}(u\_{tt})\_j^n(t-t\_n)^2, \quad (\mathbf{x},t) \in (\mathbf{SE})\_j^n \end{aligned} \tag{5.2}$$

$$f(\mathbf{x}, t) = f\_j^{\boldsymbol{n}} + (f\_\mathbf{x})\_j^{\boldsymbol{n}} (\mathbf{x} - \mathbf{x}\_j) + (f\_t)\_j^{\boldsymbol{n}} (t - t\_\mathbf{n}) + (f\_{xt})\_j^{\boldsymbol{n}} (\mathbf{x} - \mathbf{x}\_j) (t - t\_\mathbf{n})$$

$$+ \frac{1}{2} (f\_{xx})\_j^{\boldsymbol{n}} (\mathbf{x} - \mathbf{x}\_j)^2 + \frac{1}{2} (f\_{tt})\_j^{\boldsymbol{n}} (t - t\_\mathbf{n})^2. \quad (\mathbf{x}, t) \in (\mathbf{SE})\_j^{\boldsymbol{n}} \tag{5.3}$$

where *uxt*, *uxx*, *utt*, *f xt*, *f xx*, and *f tt* are the second-order derivatives of *u* and *f* , respectively. These derivatives are assumed constant within the SE. With the aid of the Eqs. (5.1) and (5.2), then one has

$$(u\_t)\_j^n = -a(u\_x)\_j^n, \quad (u\_{xt})\_j^n = -a(u\_{xx})\_j^n, \quad (u\_{tt})\_j^n = a^2(u\_{xx})\_j^n. \tag{5.4}$$

It implies that *u*, *ux*, *uxx* are the only independent variables associated with each grid point. Substituting Eqs. (5.2) and (5.3) into the integral form of Eq. (5.1). The following algebraic relation can be derived as

$$(u\_j^n \Delta x + (u\_{xx})\_j^n \frac{\Delta x^3}{24} = U\_L + U\_R + F\_L - F\_R,\tag{5.5}$$

where

#### 5.1 Construction of High-Order CESE Schemes 59

$$\begin{array}{l} U\_{L} = \frac{\Delta x}{2} u\_{j-1/2}^{n-1/2} + \frac{\Delta x^{2}}{8} (u\_{x})\_{j-1/2}^{n-1/2} + \frac{\Delta x^{3}}{48} (u\_{xx})\_{j-1/2}^{n-1/2} \\ U\_{R} = \frac{\Delta x}{2} u\_{j+1/2}^{n-1/2} - \frac{\Delta x^{2}}{8} (u\_{x})\_{j+1/2}^{n-1/2} + \frac{\Delta x^{3}}{48} (u\_{xx})\_{j+1/2}^{n-1/2} \\ F\_{L} = \frac{\Delta t}{2} f\_{j-1/2}^{n-1/2} + \frac{\Delta t^{2}}{8} (f\_{t})\_{j-1/2}^{n-1/2} + \frac{\Delta t^{3}}{48} (f\_{tt})\_{j-1/2}^{n-1/2} \\ F\_{R} = \frac{\Delta t}{2} f\_{j+1/2}^{n-1/2} + \frac{\Delta t^{2}}{8} (f\_{t})\_{j+1/2}^{n-1/2} + \frac{\Delta t^{3}}{48} (f\_{tt})\_{j+1/2}^{n-1/2} \end{array} \tag{5.6}$$

In Eqs. (5.5) and (5.6), the values at *t* = *tn*−1/2 are already known. It is necessary to calculate (*uxx* ) *n j* at *t* =*tn* before explicitly computing *u<sup>n</sup> <sup>j</sup>*. The second-order derivative, (*uxx* ) *n <sup>j</sup>*, is approximated by a central difference as

$$(\mu\_{xx})\_j^n = \frac{(\mu\_x)\_{j+1/2}^n - (\mu\_x)\_{j-1/2}^n}{\Delta x}.\tag{5.7}$$

Then *u<sup>n</sup> <sup>j</sup>*can be directly computed from Eq. (5.5), and (*ux* ) *n <sup>j</sup>*can be calculated in the same manner as the second-order *a-*α schemes using the weighted average function.

#### *5.1.2 Construction of a Third-Order 2D CESE Scheme on Uniform Mesh*

Recall the 2D scalar hyperbolic conservation law

$$
\frac{\partial \mu}{\partial t} + \frac{\partial f}{\partial x} + \frac{\partial g}{\partial y} = 0. \tag{5.8}
$$

which can be cast into the integral form

$$\oint\_{S(\text{CE}(P'))} \mathbf{h} \cdot \mathbf{n} \, \text{d}S = 0 \tag{5.9}$$

where *h* = ( *f*, *g*, *u*)is the space–time flux vector, *n* is the unit outward normal vector on the surface of the control volume. Inspired by the construction of the high-order 1D scheme, second-order Taylor expansion is utilized inside the SE to estimate *u*, *f* , and *g*. The definitions of SE and CE are kept the same as the second-order schemes (Fig. 4.1). For a point (*x*, *y*, *t*) that belongs to SE(*P'*),

$$
\mu(\mathbf{x}, \mathbf{y}, t) = \mu(\delta \mathbf{x}, \delta \mathbf{y}, \delta t)\_{P'}, \quad (\mathbf{x}, \mathbf{y}, t) \in \mathrm{SE}(P') \tag{5.10}
$$

$$f(\mathbf{x}, \mathbf{y}, t) = f(\delta \mathbf{x}, \delta \mathbf{y}, \delta t)\_{P'}, \quad (\mathbf{x}, \mathbf{y}, t) \in \mathrm{SE}(P') \tag{5.11}$$

$$\mathbf{g}(\mathbf{x}, \mathbf{y}, t) = \mathbf{g}(\delta \mathbf{x}, \delta \mathbf{y}, \delta t)\_{P'}, \quad (\mathbf{x}, \mathbf{y}, t) \in \mathrm{SE}(P') \tag{5.12}$$

where *X* (δ*x*, δ*y*, δ*t*)*N* denotes the second-order Taylor expansion about point *N* as

$$\begin{split} X(\delta x, \delta y, \delta t)\_N &= X\_N + (X\_x)\_N \delta x + \left(X\_y\right)\_N \delta y + (X\_t)\_N \delta t \\ &\quad + \frac{1}{2} (X\_{xx})\_N (\delta x)^2 + \frac{1}{2} (X\_{yy})\_N (\delta y)^2 + \frac{1}{2} (X\_{tt})\_N (\delta t)^2 \\ &\quad + \left(X\_{xy}\right)\_N \delta x \delta y + \left(X\_{yt}\right)\_N \delta y \delta t + (X\_{xt})\_N \delta x \delta t. \end{split} \tag{5.13}$$

and

$$
\delta \mathbf{x} = \mathbf{x} - \mathbf{x}\_N, \ \delta \mathbf{y} = \mathbf{y} - \mathbf{y}\_N, \ \delta t = t - t\_N. \tag{5.14}
$$

Substituting Eqs. (5.10)–(5.12) into (5.8), one obtains

$$\begin{aligned} u\_t &= -f\_x - \mathbf{g}\_y \\ u\_{tt} &= -f\_{xt} - \mathbf{g}\_{yt} \\ u\_{xt} &= -f\_{xx} - \mathbf{g}\_{xy} \\ u\_{yt} &= -f\_{xy} - \mathbf{g}\_{yy} \end{aligned} \tag{5.15}$$

Together with the chain rule, to compute the derivatives in Eq. (5.13), the only unknowns are *u*, *ux* , *uy* , *ux x* , *uyy* , and *ux y* . By integrating over the surface of the CE(*P'*), Eq. (5.9) leads to

$$\begin{aligned} \oint\_{S(V)} \mathbf{h} \cdot \mathbf{n} \, \mathrm{d}S &= \int\_{A'B'C'D'} \mathbf{h} \cdot \mathbf{n} \, \mathrm{d}S + \int\_{ABCD} \mathbf{h} \cdot \mathbf{n} \, \mathrm{d}S + \int\_{ABB'A'} \mathbf{h} \cdot \mathbf{n} \, \mathrm{d}S \\ &+ \int\_{BCC'B'} \mathbf{h} \cdot \mathbf{n} \, \mathrm{d}S + \int\_{CDB'C'} \mathbf{h} \cdot \mathbf{n} \, \mathrm{d}S + \int\_{DAA'D'} \mathbf{h} \cdot \mathbf{n} \, \mathrm{d}S = 0. \end{aligned} \tag{5.16}$$

Thereafter, with the aid of Eqs. (5.10) –(5.12) and (5.16) can be rearranged as

$$(u)\_{p'} + \frac{\left(\Delta x\right)^2}{24} (u\_{xx})\_{p'} + \frac{\left(\Delta y\right)^2}{24} (u\_{yy})\_{p'} = \frac{1}{4} \left(\overline{u} + \frac{\Delta t}{\Delta x} \overline{f} + \frac{\Delta t}{\Delta y} \overline{g}\right),\qquad(5.17)$$

where

$$\begin{split} \bar{u} &= \hat{u}\left(\frac{\Delta x}{4}, \frac{\Delta \mathbf{y}}{4}, 0\right)\_A + \hat{u}\left(-\frac{\Delta x}{4}, \frac{\Delta \mathbf{y}}{4}, 0\right)\_B \\ &+ \hat{u}\left(-\frac{\Delta x}{4}, -\frac{\Delta \mathbf{y}}{4}, 0\right)\_C + \hat{u}\left(\frac{\Delta x}{4}, -\frac{\Delta \mathbf{y}}{4}, 0\right)\_D \end{split} \tag{5.18}$$
 
$$\bar{f} = \hat{f}\left(0, \frac{\Delta \mathbf{y}}{4}, \frac{\Delta t}{4}\right)\_A - \hat{f}\left(0, \frac{\Delta \mathbf{y}}{4}, \frac{\Delta t}{4}\right)\_B$$

$$-\left.\hat{f}\left(0, -\frac{\Delta \mathbf{y}}{4}, \frac{\Delta t}{4}\right)\_C + \hat{f}\left(0, -\frac{\Delta \mathbf{y}}{4}, \frac{\Delta t}{4}\right)\_D\right. \tag{5.19}$$

$$
\begin{split}
\bar{g} &= \hat{g}\left(\frac{\Delta x}{4}, 0, \frac{\Delta t}{4}\right)\_A + \hat{g}\left(-\frac{\Delta x}{4}, 0, \frac{\Delta t}{4}\right)\_B \\ &- \hat{g}\left(-\frac{\Delta x}{4}, 0, \frac{\Delta t}{4}\right)\_C - \hat{g}\left(\frac{\Delta x}{4}, 0, \frac{\Delta t}{4}\right)\_D
\end{split}
\tag{5.20}
$$

and

$$\begin{split} \hat{u}(\delta x, \delta \mathbf{y}, \delta t)\_N &= (u)\_N + (u\_x)\_N \delta x + \left(u\_\mathbf{y}\right)\_N \delta \mathbf{y} + (u\_t)\_N \delta t \\ &+ \frac{1}{6} (u\_{xx})\_N (\delta x)^2 + \frac{1}{6} (u\_{\mathbf{y}\mathbf{y}})\_N (\delta \mathbf{y})^2 + \left(u\_{\mathbf{x}\mathbf{y}}\right)\_N \delta x \delta \mathbf{y}, \end{split} \tag{5.21}$$

$$\begin{split} \hat{f}(\delta \mathbf{x}, \delta \mathbf{y}, \delta t)\_N &= (f)\_N + (f\_\mathbf{x})\_N \delta \mathbf{x} + \left(f\_\mathbf{y}\right)\_N \delta \mathbf{y} + (f\_\mathbf{t})\_N \delta t \\ &+ \frac{1}{6} (f\_\mathbf{y})\_N (\delta \mathbf{y})^2 + \frac{1}{6} (f\_\mathbf{t})\_N (\delta t)^2 + \left(f\_\mathbf{y}\right)\_N \delta y \delta t, \end{split} \tag{5.22}$$

$$\begin{split} \hat{\mathbf{g}}(\delta \mathbf{x}, \delta \mathbf{y}, \delta t)\_N &= (\mathbf{g})\_N + (\mathbf{g}\_x)\_N \delta \mathbf{x} + \left( \mathbf{g}\_y \right)\_N \delta \mathbf{y} + (\mathbf{g}\_t)\_N \delta t \\ &+ \frac{1}{6} (\mathbf{g}\_{xx})\_N (\delta \mathbf{x})^2 + \frac{1}{6} (\mathbf{g}\_{tt})\_N (\delta t)^2 + (\mathbf{g}\_{xt})\_N \delta \mathbf{x} \delta t. \end{split} \tag{5.23}$$

Prior to explicitly computing (*u*)*p* from Eq. (5.17), (*uxx* )*p* and *uyy <sup>p</sup>* are required to be obtained from the finite difference of the interpolations from the previous time step. The second-order derivative (*uxx* )*p* can be estimated as

$$\left(u\_{xx}\right)\_{p'}^+ = \frac{(u\_x)\_B + \frac{\Delta t}{2}(u\_{xt})\_B - (u\_x)\_A - \frac{\Delta t}{2}(u\_{xt})\_A}{\Delta x},\tag{5.24}$$

$$(\mu\_{xx})\_{p'}^- = \frac{(\mu\_x)\_C + \frac{\Delta t}{2}(\mu\_{xt})\_C - (\mu\_x)\_D - \frac{\Delta t}{2}(\mu\_{xt})\_D}{\Delta x}. \tag{5.25}$$

(*uxx* )*p* is then computed from simple average as

$$(u\_{xx})\_{p'} = \frac{(u\_{xx})\_{p'}^+ + (u\_{xx})\_{p'}^-}{2}. \tag{5.26}$$

Similarly, the derivative *uyy <sup>p</sup>* is computed as

$$\left(\left(u\_{\rm yr}\right)\right)\_{p'}^+ = \frac{\left(u\_{\rm y}\right)\_D + \frac{\Delta t}{2}\left(u\_{\rm yt}\right)\_D - \left(u\_{\rm y}\right)\_A - \frac{\Delta t}{2}\left(u\_{\rm yt}\right)\_A}{\Delta \mathbf{y}},\tag{5.27}$$

$$\left(\left(\mu\_{\rm yy}\right)\_{p'}^{-} = \frac{\left(\mu\_{\rm y}\right)\_{\rm C} + \frac{\Delta t}{2}\left(\mu\_{\rm yt}\right)\_{\rm C} - \left(\mu\_{\rm y}\right)\_{\rm B} - \frac{\Delta t}{2}\left(\mu\_{\rm yt}\right)\_{\rm B}}{\Delta \rm y},\tag{5.28}$$

$$\left(\mu\_{\rm yr}\right)\_{p'} = \frac{\left(\mu\_{\rm yr}\right)\_{p'}^+ + \left(\mu\_{\rm yr}\right)\_{p'}^-}{2}.\tag{5.29}$$

Equations (5.17), (5.26), and (5.29) form the procedure to get the values of (*u*)*p* , (*uxx* )*p* and *uyy <sup>p</sup>* at the new time level. To complete the time marching, we need to further compute *ux* , *uy* , and *ux y* . The weighted average of the sided-estimations computes the first-order derivatives as

$$\left(\left(u\_x\right)\_{p'}^+ = \frac{1}{\Delta x} \left[ u\left(0, 0, \frac{\Delta t}{2}\right)\_B + u\left(0, 0, \frac{\Delta t}{2}\right)\_C - 2\left(u\_x\right)\_{p'} \right],\tag{5.30}$$

$$(u\_x)\_{p'}^- = \frac{1}{\Delta x} \left[ 2(u\_x)\_{p'} - \mu \left( 0, 0, \frac{\Delta t}{2} \right)\_A - \mu \left( 0, 0, \frac{\Delta t}{2} \right)\_D \right],\tag{5.31}$$

$$\left(\left(u\_{\rm y}\right)\_{p'}^+ = \frac{1}{\Delta y} \Big[ u\left(0, 0, \frac{\Delta t}{2}\right)\_C + u\left(0, 0, \frac{\Delta t}{2}\right)\_D - 2\left(u\_{\rm x}\right)\_{p'} \Big],\tag{5.32}$$

$$\left(\left(\mu\_{\rm y}\right)\_{p'}^{-} = \frac{1}{\Delta y} \left[ 2\left(\mu\_{\rm y}\right)\_{p'} - \mu\left(0, 0, \frac{\Delta t}{2}\right)\_{A} - \mu\left(0, 0, \frac{\Delta t}{2}\right)\_{B} \right],\tag{5.33}$$

$$(\mu\_x)\_{p'} = W(\left(\mu\_x^-\right)\_{p'}, \left(\mu\_x^+\right)\_{p'}, \alpha), \tag{5.34}$$

$$\left(\left(\mu\_{\rm y}\right)\_{P'} = W\left(\left(\mu\_{\rm y}^{-}\right)\_{P'}, \left(\mu\_{\rm y}^{+}\right)\_{P'}, \alpha\right). \tag{5.35}$$

Meanwhile, the mixed derivative *uxy* is computed from

$$\left(\left(u\_{xy}\right)\_{p'}^+ = \frac{(u\_x)\_D + \frac{\Delta t}{2}(u\_{xt})\_D - (u\_x)\_A - \frac{\Delta t}{2}(u\_{xt})\_A}{\Delta \mathbf{y}},\tag{5.36}$$

$$\left(\mu\_{xy}\right)^{-}\_{p'} = \frac{(\mu\_x)\_C + \frac{\Delta t}{2}(\mu\_{xt})\_C - (\mu\_x)\_B - \frac{\Delta t}{2}(\mu\_{xt})\_B}{\Delta \mathbf{y}},\tag{5.37}$$

$$\left(\mu\_{\rm yx}\right)\_{p'}^+ = \frac{(\mu\_x)\_B + \frac{\Delta t}{2}(\mu\_{xt})\_B - (\mu\_x)\_A - \frac{\Delta t}{2}(\mu\_{xt})\_A}{\Delta x},\tag{5.38}$$

$$\left(\mu\_{\rm yx}\right)^{-}\_{p'} = \frac{(\mu\_x)\_C + \frac{\Delta t}{2}(\mu\_{xt})\_C - (\mu\_x)\_D - \frac{\Delta t}{2}(\mu\_{xt})\_D}{\Delta x},\tag{5.39}$$

$$\left(\left(\mu\_{\rm xy}\right)\_{p'} = \frac{1}{4}\left[\left(\mu\_{\rm xy}\right)\_{p'}^{+} + \left(\mu\_{\rm xy}\right)\_{p'}^{-} + \left(\mu\_{\rm yx}\right)\_{p'}^{+} + \left(\mu\_{\rm yx}\right)\_{p'}^{-}\right].\tag{5.40}$$

#### *5.1.3 Construction of a Third-Order 2D CESE Scheme on Unstructured Mesh*

Similar to constructing the third-order scheme on a a uniform mesh, second-order Taylor expansion is utilized here, and the variables stored at each grid point are still *u*, *ux* , *uy* , *ux x* , *uyy* , and *ux y* . Now we recollect the definitions of SE and CE in the 2D CESE scheme on the unstructured mesh as depicted in Fig. 4.4, and the flux balancing relation was derived in Eq. (4.48). From the summation of the conservation law for all the sub-CEs, we arrive at the flux balancing relation for CE(*V i* ):

$$\iint\limits\_{\mathbf{C}\_{1}^{\prime}O\_{1}^{\prime}\cdots C\_{M}O\_{M}^{\prime}} u \, d\sigma = \sum\_{m=1}^{M} \iint\limits\_{\mathbf{C}\_{m}O\_{m}^{\prime}V\_{i}O\_{m-1}} u \, d\sigma - \sum\_{m=1}^{M} \iint\limits\_{\mathbf{C}\_{m}O\_{m}^{\prime}O\_{m}^{\prime}C\_{m}^{\prime} + O\_{m-1}C\_{m}C\_{m}^{\prime}O\_{m-1}^{\prime}} \mathbf{V} \cdot \mathbf{n} \, d\sigma. \tag{5.41}$$

The LHS of Eq. (5.41) represents the integration of *u* on the top surface of CE(*V i* ), i.e.,

$$\begin{split} \iint \limits\_{\mathbf{C}'\_1 \mathbf{C}'\_1 \dots \mathbf{C}'\_M} u \, d\sigma &= u \big( \mathbf{G}'\_l \big) \cdot \mathbf{S} + \text{TERM1} = u \big( \mathbf{G}'\_l \big) \cdot \mathbf{S} + \sum\_{m=1}^M \iint \limits\_{\mathbf{C}'\_m \mathbf{C}'\_m \mathbf{G}'\_l \big( \mathbf{G}'\_{m-1} \big)} \frac{1}{2} \frac{\partial^2 u \big( \mathbf{G}'\_l \big)}{\partial x^2} (\delta x)^2 \\ &+ \frac{1}{2} \frac{\partial^2 u \big( \mathbf{G}'\_l \big)}{\partial y^2} (\delta y)^2 + \frac{\partial^2 u \big( \mathbf{G}'\_l \big)}{\partial x \partial y} (\delta x) (\delta y) dx dy, \end{split} \tag{5.42}$$

where δ*x* = *x* − *x*(*G <sup>i</sup>*), δ*y* = *y* − *y*(*G <sup>i</sup>*), and *S* represents the area of the polygon *C* 1*O* 1 ...*C M O <sup>M</sup>* . The second term in RHS (TERM1) of Eq. (5.42) denotes the integration of the high-order terms of *u*. The two terms in RHS of Eq. (5.41) can be respectively expressed as

$$\sum\_{m=-1}^{M} \iint\_{\mathbf{C}\_{\text{m}} \otimes \mathbf{C}\_{\text{l}} \otimes \mathbf{C}\_{\text{m}-1}} u \, d\sigma = \text{TERM2} = \sum\_{m=1}^{M} \iint\_{\mathbf{C}\_{\text{m}} \otimes \mathbf{C}\_{\text{l}} \otimes \mathbf{C}\_{\text{m}-1}} u(\mathbf{C}\_{\text{m}}) + \frac{\partial u(\mathbf{C}\_{\text{m}})}{\partial \mathbf{x}} (\delta \mathbf{x}) + \frac{\partial u(\mathbf{C}\_{\text{m}})}{\partial \mathbf{y}} (\delta \mathbf{y})$$

$$+ \frac{1}{2} \frac{\partial^{2} u(\mathbf{C}\_{\text{m}})}{\partial x^{2}} (\delta x)^{2} + \frac{1}{2} \frac{\partial^{2} u(\mathbf{C}\_{\text{m}})}{\partial y^{2}} (\delta \mathbf{y})^{2} + \frac{\partial^{2} u(\mathbf{C}\_{\text{m}})}{\partial x \partial \mathbf{y}} (\delta x) (\delta \mathbf{y}) d\mathbf{x} d\mathbf{y},\tag{5.43}$$

where δ*x* = *x* − *x*(*Cm*), δ*y* = *y* − *y*(*Cm*).

$$\sum\_{m=1}^{M} \iint\_{\mathbf{C}\_{m}\mathbf{O}\_{m}\mathbf{O}\_{m}^{\prime}\mathbf{C}\_{m}^{\prime} + \mathbf{O}\_{n-1}\mathbf{C}\_{n}\mathbf{C}\_{m}^{\prime}\mathbf{O}\_{n-1}^{\prime}} \mathbf{V} \cdot \mathbf{n} d\sigma = \text{TERM3} = \sum\_{m=1}^{M} \left[ FLUX\_{1}^{(m)} + FLUX\_{2}^{(m)} \right]. \tag{5.44}$$

with

*FLU X*(*m*) <sup>1</sup><sup>=</sup> *<sup>t</sup>* <sup>2</sup> *<sup>y</sup>* (*m*) 1 <sup>2</sup> *a*=0 2 −*a b*=0 2− *a*−*<sup>b</sup> c*=0 1 (*a* + *b* + 1)*a*!*b*!(*c* + 1)! ∂*a*+*b*+*c f* (*Cm* ) ∂*xa* ∂*yb*∂*tc x*(*m*) 1 *<sup>a</sup> y* (*m*) 1 *b t*  2 *c*  <sup>−</sup> *<sup>t</sup>* <sup>2</sup> *x*(*m*) 1 <sup>2</sup> *a*=0 2 −*a b*=0 2− *a*−*<sup>b</sup> c*=0 1 (*a* + *b* + 1)*a*!*b*!(*c* + 1)! ∂*a*+*b*+*cg*(*Cm* ) ∂*xa* ∂*yb*∂*tc x*(*m*) 1 *<sup>a</sup> y* (*m*) 1 *b t*  2 *c*  , (5.45) *FLU X*(*m*) <sup>2</sup><sup>=</sup> *<sup>t</sup>* <sup>2</sup> *<sup>y</sup>* (*m*) 2 <sup>2</sup> *a*=0 2 −*a b*=0 2− *a*−*<sup>b</sup> c*=0 1 (*a* + *b* + 1)*a*!*b*!(*c* + 1)! ∂*a*+*b*+*c f* (*Cm* ) ∂*xa* ∂*yb*∂*tc* −*x*(*m*) 2 *<sup>a</sup>* −*y* (*m*) 2 *b t*  2 *c*  <sup>−</sup> *<sup>t</sup>* <sup>2</sup> *x*(*m*) 2 <sup>2</sup> *a*=0 2 −*a b*=0 2− *a*−*<sup>b</sup> c*=0 1 (*a* + *b* + 1)*a*!*b*!(*c* + 1)! ∂*a*+*b*+*cg*(*Cm* ) ∂*xa* ∂*yb*∂*tc* −*x*(*m*) 2 *<sup>a</sup>* −*y* (*m*) 2 *b t*  2 *c*  , (5.46)

in which *x*(*m*) <sup>1</sup><sup>=</sup> *<sup>x</sup>*(*Om*) <sup>−</sup> *<sup>x</sup>*(*Cm*), *y*(*m*) <sup>1</sup><sup>=</sup> *<sup>y</sup>*(*Om*) <sup>−</sup> *<sup>y</sup>*(*Cm*), *x*(*m*) <sup>2</sup>= *x*(*Cm*) − *<sup>x</sup>*(*Om*−<sup>1</sup>), *y*(*m*) <sup>2</sup>= *y*(*Cm*)− *y*(*Om*−<sup>1</sup>). *FLUX(m)* 1 and *FLUX(m)* 2 represent the fluxes across the adjacent CE surfaces on which the points are defined within SE(*Cm*).

By substituting Eqs. (5.42)–(5.44) into (5.41) and moving the high-order terms to RHS. Consequently, the value of *u*(*G <sup>i</sup>*) can be expressed by

$$
\mu\left(G\_i'\right) = \frac{1}{S}(\text{TERM2} - \text{TERM3} - \text{TERM1}).\tag{5.47}
$$

Two remarks should be noted regarding Eq. (5.47). First of all, since *u*, *f* , and *g*  are nonlinearly distributed on the surfaces of CE, the integral terms appeared in Eqs. (5.42) and (5.43) are not simple multiplications of the value at the centroid and the surface area as what has been implemented in the second-order scheme. In contrast, they are now calculated using the coordinate transformation method:

$$\chi = \sum\_{i=1}^{4} \chi\_i N\_i(\xi, \ \eta),\tag{5.48}$$

$$\mathbf{y} = \sum\_{i=1}^{4} \mathbf{y}\_i N\_i(\xi, \ \eta),\tag{5.49}$$

where

$$\begin{aligned} N\_1 &= \frac{1}{4}(1-\xi)(1-\eta), & N\_2 &= \frac{1}{4}(1+\xi)(1-\eta),\\ N\_3 &= \frac{1}{4}(1+\xi)(1+\eta), & N\_4 &= \frac{1}{4}(1-\xi)(1+\eta). \end{aligned} \tag{5.50}$$

Thereafter, the integration of an arbitrary function φ(*x*, *y*) over the irregular quadrangle *A* is transformed into the integration over a rectangle (Fig. 5.1) as

**Fig. 5.1** Sketch of the coordinate transformation

$$\iint\_{A} \phi(\mathbf{x}, \mathbf{y}) ds = \int\_{-1}^{1} d\eta \int\_{-1}^{1} \phi[\mathbf{x}(\xi, \eta), \mathbf{y}(\xi, \eta)] \mathbf{J} d\xi,\tag{5.51}$$

where **J** is the Jacobian matrix of the coordinate transformation,

$$\mathbf{J} = \frac{\partial(\mathbf{x}, \mathbf{y})}{\partial(\boldsymbol{\xi}, \boldsymbol{\eta})} = \begin{vmatrix} \sum\_{i=1}^{4} \boldsymbol{\chi}\_{i} \frac{\partial \boldsymbol{N}\_{i}}{\partial \boldsymbol{\xi}} \sum\_{i=1}^{4} \boldsymbol{\chi}\_{i} \frac{\partial \boldsymbol{N}\_{i}}{\partial \boldsymbol{\xi}}\\ \sum\_{i=1}^{4} \boldsymbol{\chi}\_{i} \frac{\partial \boldsymbol{N}\_{i}}{\partial \boldsymbol{\eta}} \sum\_{i=1}^{4} \boldsymbol{\chi}\_{i} \frac{\partial \boldsymbol{N}\_{i}}{\partial \boldsymbol{\eta}} \end{vmatrix}.\tag{5.52}$$

In addition, TERM2 and TERM3 in Eq. (5.47) can be computed from the information at *t* = *tn*−1/2. TERM1 consists of high-order terms (*ux x* , *uyy* , and *ux y* ) of *u* at *t* = *tn*. As a result, the high-order derivatives must be solved prior to computing *u*. The idea to calculate the high-order derivatives are analogue to that in Sects. 5.1.1 and 5.1.2. For example, from the Taylor expansion *C <sup>m</sup>*we have

$$
\mu\_X \left( \mathcal{C}\_{m-1}^{\prime} \right) - \mu\_X \left( \mathcal{C}\_m^{\prime} \right) = \mu\_{XX} \left( \mathcal{C}\_m^{\prime} \right) \left[ x(\mathcal{C}\_{m-1}^{\prime}) - x(\mathcal{C}\_m^{\prime}) \right] + \mu\_{X\uparrow} \left( \mathcal{C}\_m^{\prime} \right) \left[ \mathbf{y}(\mathcal{C}\_{m-1}^{\prime}) - \mathbf{y}(\mathcal{C}\_m^{\prime}) \right], \tag{5.55}
$$

$$
\mu\_X \left( \mathbf{C}'\_{m+1} \right) - \mu\_X \left( \mathbf{C}'\_m \right) = \mu\_{X\mathbf{X}} \left( \mathbf{C}'\_m \right) \left[ \mathbf{x}(\mathbf{C}'\_{m+1}) - \mathbf{x}(\mathbf{C}'\_m) \right] + \mu\_{X\mathbf{Y}} \left( \mathbf{C}'\_m \right) \left[ \mathbf{y}(\mathbf{C}'\_{m+1}) - \mathbf{y}(\mathbf{C}'\_m) \right]. \tag{5.54}
$$

On the other hand, *ux* (*C <sup>m</sup>*−<sup>1</sup>), *ux* (*C m*), *ux* (*C <sup>m</sup>*+<sup>1</sup>) can be interpolated from the previous time level,

$$
\mu\_x \left( C\_{m-1}' \right) = \mu\_x \left( C\_{m-1} \right) + \mu\_{xx} \left( C\_{m-1} \right) \frac{\Delta t}{2}
$$

$$
\mu\_x \left( C\_m' \right) = \mu\_x \left( C\_m \right) + \mu\_{xx} \left( C\_m \right) \frac{\Delta t}{2}
$$

$$
\mu\_x \left( C\_{m+1}' \right) = \mu\_x \left( C\_{m+1} \right) + \mu\_{xx} \left( C\_{m+1} \right) \frac{\Delta t}{2} \tag{5.55}
$$

By substituting Eq. (5.55) into Eqs. (5.53) and (5.54), a pair of equations with two unknowns (*uxx* (*C <sup>m</sup>*), *ux y* (*C m*)) are ready to be solved by Cramer's rule. A simple average (or weighted average to suppress oscillations near discontinuities) can be applied for *M* pairs of derivatives. Once all the second-order derivatives are obtained under the same routine, *u*(*G <sup>i</sup>*) can be immediately calculated from Eq. (5.47). The last unknowns remain to be computed are *ux* and *uy* , which can be calculated similarly in Sect. 4.2, e.g.,

$$\begin{cases} u\left(G\_i'\right) + u\_x\left(G\_i'\right)\delta x + u\_y\left(G\_i'\right)\delta y + \frac{1}{2}u\_{xx}\left(G\_i'\right)\left(\delta x\right)^2 + \frac{1}{2}u\_{yy}\left(G\_i'\right)\left(\delta y\right)^2\\ + u\_{xy}\left(G\_i'\right)\left(\delta x\right)\left(\delta y\right) = u\left(C\_m\right) + u\_t\left(C\_m\right)\frac{\Delta t}{2} + \frac{1}{2}u\_{tt}\left(C\_m\right)\left(\frac{\Delta t}{2}\right)^2 \end{cases},\quad(5.56)$$

where δ*x* = *x*(*Cm*)−*x*(*G <sup>i</sup>*), δ*y* = *y*(*Cm*)− *y*(*G <sup>i</sup>*). Similarly, on applying Eq. (5.56) to *Cm*+1, a pair of *ux* and *uy* can be easily determined. Then, we apply the weighted average function to obtain the optimal derivatives. Finally, *u*(*V <sup>i</sup>*) and its derivatives are interpolated by Taylor expansion from the information at *u*(*G i* ).

#### **5.2 Numerical Examples**

The Shu-Osher problem consists of a Mach 3 shock interacting with an entropy wave, which requires the capability of capturing fine structures in the flow. The computational domain [0, 10] is discretized by a uniform grid size of 0.05, the nonreflection boundary condition is imposed on both sides, and the initial condition is described as

$$\mathbf{x}(\rho,\ u,\ p) = \begin{cases} (3.857143, \ 2.629369, \ 10.33333) \ x \le 1 \\ (1 + 0.2\sin(5x), \ 0, \ 1) \end{cases} \tag{5.57}$$

The density distributions at *t* = *1.8* are depicted in (Fig. 5.2). Compared to the benchmark solution, which is calculated with *a*-α scheme using very high grid resolution, higher-order CESE schemes are less dissipative and can reveal better flow structures.

The second example is the Mach 3 wind tunnel problem with step. This problem has been described in Sect. 4.3. Here, we use the identical numerical setups except with a mesh size of 1/250. As depicted in Fig. 5.3, vortex-like structures generated in the results using high-order schemes, however, disappeared in the low-order schemes due to the increasing dissipation. It is not surprising that the triangular mesh result surpasses the quadrilateral ones because the mesh area is only half of the latter when the mesh sizes are the same.

**Fig. 5.2** Density distributions for Shu-Osher problem at *t* = 1.8. **a** Entire view, **b** Enlargement. Courtesy of Shen [7]

**Fig. 5.3** Density contours for Mach 3 wind tunnel with a step at *t* = 4.0. Courtesy of Shen [7]

#### **References**


In *48th AIAA Aerospace Sciences Meeting Including the New Horizons Forum and Aerospace Exposition.* 


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **Chapter 6 Numerical Features of CESE Schemes**

In this chapter, some remarks are made on the numerical characteristics of the CESE schemes described in foregoing chapters. Due to the special formulation, rigorous analysis of the CESE schemes will take more efforts than that of traditional finite difference schemes. However, it is possible to extend the widely used modified equation analysis, modified wavenumber analysis, and von Neumann stability analysis to the CESE schemes. Here, a stability analysis for the second-order upwind CESE scheme applied to the linear scalar convection equation is presented. Additionally, a rough estimation of the algorithm complexity is presented to show the computational efficiency of the CESE method. In order to demonstrate the accuracy of the CESE method in comparison with established numerical approaches, simulation results of some canonical problems are shown.

#### **6.1 Stability**

All CESE schemes in this book are explicit time-stepping schemes. Similar to other explicit schemes, the numerical stability of CESE schemes is closely related to the Courant–Friedrichs–Lewy (CFL) condition. Only when the CFL condition is satisfied, the system can remain numerically stable. Consequently, the time-step size used by a CESE scheme has to be restricted according to the upper bound of the stable CFL number ν for that scheme.

It has been proven in [1, 2] that the stability conditions for the second-order *a*  scheme (Sect. 2.3) and the second-order *a*–α scheme (Sect. 3.1) are both ν ≤ 1. For the second-order upwind CESE scheme (Sect. 3.3), the stable upper bound of the CFL number is also unity. This stability analysis is first given by Jiang et al. [3] and introduced as follows.

For illustrative purpose, we consider the simplest upwind CESE scheme applied to the 1D linear scalar convection equation

70 6 Numerical Features of CESE Schemes

$$\frac{\partial u}{\partial t} + a \frac{\partial u}{\partial x} = 0 \tag{6.1}$$

where *a* is a positive constant. Without using any slope limiter, the upwind CESE scheme for Eq. (6.1) can be written in a matrix form as

$$\mathbf{q}\_{j}^{n} = \mathbf{Q}\_{L}\mathbf{q}\_{j-1/2}^{n-1/2} + \mathbf{Q}\_{R}\mathbf{q}\_{j+1/2}^{n-1/2} \tag{6.2}$$

with the vector of time-marching variables

$$\mathbf{q}\_{j}^{n} = \begin{bmatrix} \boldsymbol{u}\_{j}^{n} \\ \boldsymbol{\Delta x} \\ \hline \boldsymbol{4} \end{bmatrix}, \tag{6.3}$$

and the coefficient matrices (recall Chap. 3)

$$\mathcal{Q}\_L = \frac{1}{2} \begin{bmatrix} 1+\nu & 1-\nu^2 \\ -1+\nu & -1+4\nu-\nu^2 \end{bmatrix}, \qquad \mathcal{Q}\_R = \frac{1}{2} \begin{bmatrix} 1-\nu & -1+\nu^2 \\ 1-\nu & -1+\nu^2 \end{bmatrix}, \tag{6.4}$$

where ν is the CFL number, which is defined as

$$
\upsilon \equiv a \frac{\Delta t}{\Delta x}. \tag{6.5}
$$

Following the routine of von Neumann stability analysis, we assume the vector *q*  to be a Fourier component with arbitrary wavenumber *k* such that

$$\mathbf{q}\_{j}^{n} = A\_{n}(\theta)e^{ij\theta}, \ (-\pi \prec \theta \le \pi), \tag{6.6}$$

where θ = *k*Δ*x* is the reduced wavenumber, and *i* is the imaginary unit. Substituting Eqs. (6.6) into (6.2) yields

$$A\_n(\theta) = \left[ \mathcal{Q}\_L e^{-i\theta/2} + \mathcal{Q}\_R e^{i\theta/2} \right] \mathbf{A}\_{n-1/2}(\theta). \tag{6.7}$$

Hence, the amplification matrix for each half step is

$$\begin{split} \mathbf{M} &= \mathbf{Q}\_L e^{-i\theta/2} + \mathbf{Q}\_R e^{i\theta/2} \\ &= \begin{bmatrix} \cos\frac{\theta}{2} - i\nu\sin\frac{\theta}{2} & i(\nu^2 - 1)\sin\frac{\theta}{2} \\ i(1 - \nu)\sin\frac{\theta}{2} & (2\nu - 1)\cos\frac{\theta}{2} + i(\nu^2 - 2\nu)\sin\frac{\theta}{2} \end{bmatrix} .\end{split} \tag{6.8}$$

The spectral radius of matrix *M*, denoted as ρ(ν, θ), is calculated numerically as a function of the CFL number (0 ≤ ν ≤ 1) and reduced wavenumber (−π ≤ θ ≤ π). In Fig. 6.1, it can be seen that ρ(ν, θ) ≤ 1 holds for all wavenumbers as long as 0 ≤ ν ≤ 1. Thus, this second-order upwind CESE scheme is numerically stable when ν

**Fig. 6.1** Spectral radius ρ(ν, θ) of matrix *M* in the von Neumann analysis of the upwind CESE scheme. Courtesy of Jiang [3]

is less than or equal to 1. Also shown in Fig. 6.1, as the CFL number approaches 0 or 1, the numerical dissipation in the second-order upwind CESE scheme vanishes.

For high-order CESE schemes, the stability analysis is more complicate and the stable upper bound of the CFL number depends on the scheme's construction. A novel approach to construct highly stable high-order CESE schemes was proposed by Chang [4], using the same space–time stencil as that in the second-order CESE scheme. The same CFL-number constraint for numerical stability (i.e., ν ≤ 1) is retained, which is favourable for explicit time-stepping computations.

#### **6.2 Efficiency**

The CESE method has a highly compact stencil in space and time, regardless of the order of the scheme. This compactness feature provides the CESE method a great advantage when parallel computation is performed.

For a CESE scheme, the nominal order of accuracy and spatial dimensionality will determine the number of independent variables (denoted by *K*) to be stored and updated at each solution point. For a *M*-th-order accurate (*M* ≥ 2) and spatially *D*-dimensional (*D* = 1, 2, or 3) CESE scheme, *K* can be determined by a precise counting and

$$K(D,M) = \begin{cases} M, & D=1\\ M(M+1)/2, & D=2\\ M(M+1)(M+2)/6 \ D=3 \end{cases} \tag{6.9}$$

Notably, the number *K* in the CESE scheme is the same as the number of degrees of freedom for each cell in the discontinuous Galerkin (DG) scheme. Therefore, if the CESE and DG schemes are of the same order and used to solve the same problem on the same mesh, the memory requirements for a CESE scheme and a DG scheme are comparable.

The CESE method requires neither a reconstruction procedure based on a wide stencil nor a multi-stage time-integration technique. Besides, no upwind operation related to flux calculation is included in central CESE schemes. Nevertheless, the CESE method adopts the Cauchy–Kowalewski procedure to express the temporal derivatives in terms of spatial derivatives. It also takes time to update the spatial derivatives using the algorithms described in Chaps. 3, 4 and 5. The combined effect of all these factors determines the computational speed of the CESE method.

#### **6.3 Accuracy**

Generally, the accuracy of a CESE scheme depends on the order of the Taylor expansions, i.e., the degree of the approximation polynomials that are used to approximate unknowns and fluxes within each solution element. To achieve the theoretical *M*-th order of accuracy of the scheme in space and time, the (*M*−1)-th-order Taylor expansions are needed. According to practical applications, the second-order CESE schemes, including the central and upwind schemes, can usually provide a satisfactory trade-off between accuracy and efficiency.

The ability to capture shock waves and contact discontinuities with high resolution is one of the excellent features of the CESE scheme. Here, two canonical flow problems are used to demonstrate the accuracy of CESE schemes by comparing the numerical simulation results by CESE and other schemes. As the first example, Sod's shock-tube problem [5] is solved by both the CESE method and the finite volume method (FVM). The CESE code is an implementation of the second order *a*–α scheme described in Chap. 3, while the FVM counterpart incorporates the second-order MUSCL (Monotonic Upstream-centered Scheme for Conservation Laws) reconstruction, the van Leer limiter, and the HLLC Riemann solver. A uniform mesh with 200 cells and a CFL number of 6.9 are employed in both the CESE and the FVM computations. Density and internal energy distributions from *x* = − 1 to *x* = 1 at time *t* = 6.5 computed by CESE and FVM schemes are plotted in Fig. 6.2a and b, respectively, in comparison with the exact solution. As shown in Fig. 6.2, a narrower contact discontinuity is resolved by the CESE scheme than the FVM counterpart. Meanwhile the CESE results match the exact solution better than the FVM results at the tail of the rarefaction wave.

(b) CESE result vs. FVM result: internal energy

(b) Upwind CESE result vs. FVM result

**Fig. 6.3** Comparison of FVM and CESE results for Mach 3 step problem: density contours at time *t* = 4.6. Courtesy of Jiang [3]

In the second example, the problem of a Mach-3 inviscid flow through a 2D tunnel with a forward-facing step inside is solved by three schemes (*a*–α CESE, upwind CESE, and upwind FVM) of second-order accuracy in space and time. The initial conditions, boundary conditions, and geometrical parameters are the same as that in [6]. The upwind FVM flow solver is provided by Hao et al. [7]. The same computational mesh and the same Courant–Friedrichs–Lewy (CFL) number of 6.95 are employed for all the three schemes. Density fields at time *t* = 4.0 are shown in Fig. 6.3a and b. Compared with the upwind FVM scheme, both the *a*–α and upwind CESE schemes can accurately capture the structure of shock waves, but only the upwind CESE scheme can resolve the phenomenon of shear-layer instability (see the lower part of Fig. 6.3b), which is not observed in the present *a*-α CESE and FVM results.

To assess the accuracy of the high-order CESE schemes described in Chap. 5, the Shu–Osher problem and the double Mach-reflection problem were studied intensively in Ref. [8], where the fourth-order CESE scheme provided a good resolution of fine structures in the flows.

#### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **Chapter 7 Application: Compressible Multi-fluid Flows**

Multi-fluid flows involving shock-accelerated inhomogeneities and shock-induced instability play essential roles in a wide variety of problems including, but not limited to, supersonic combustion [1], inertial confinement fusion [2], and supernova explosion [3]. Numerical simulations of these complex flows prove to be challenging in the presence of moving and deformable material interfaces, especially for fluids with large differences in their densities or thermodynamic properties. Therefore, a discontinuity-capturing, mass-conserving, and positivity-preserving scheme is desirable for compressible multi-fluid simulations.

Qamar et al. [4] implemented the CESE method of Chang [5] and Zhang et al. [6] for solving the one- and two-dimensional compressible two-fluid models of Kreeft and Koren [7]. Numerical simulations were performed for gas–liquid Riemann problems and interactions of air shock waves with inhomogeneities containing lighter and heavier gases. In comparison with the non-oscillatory central scheme [5] and the kinetic flux-vector splitting scheme [8], the CESE scheme gives better resolution of discontinuities. Because of its outstanding ability for discontinuity capturing, the studies on Richtmyer–Meshkov instability (RMI) [9, 10] and Rayleigh–Taylor instability (RTI) [11, 12] have greatly benefited from the development of this method. Those numerical results not only eliminated disturbing factors in experiments (e.g., dimensional effect, boundary effect, liquid mist scattering), but also provided more ideal working conditions (that cannot be easily realized in experiment) and more flow details (that cannot be timely recorded by camera).

Recently, the upwind time–space CESE method was used extensively in studies of shock-accelerated inhomogeneous flows [13–15] and RMI [16–18]. The compressible two-fluid flows are described by a volume-fraction-based five-equation model [19] coupled with the stiffened gas equation of state [20]. Extensive numerical simulations were carried out using the maximum-principle-satisfying upwind CESE scheme [13], which is an improved version of the upwind CESE scheme presented in Chap. 4. The maximum-principle satisfying property is achieved by adopting a very simple limiter proposed by [21]. Furthermore, the ability to capture contact discontinuities

C.-Y. Wen et al., *Space–Time Conservation Element and Solution Element Method*, Engineering Applications of Computational Methods 13, https://doi.org/10.1007/978-981-99-0876-9\_7

(material interfaces) can be enhanced by employing the HLLC Riemann solver in the upwind procedure. As a result, challenging numerical simulations of the gas–gas and gas–liquid Riemann problems were successfully performed, in which the mass of each fluid component is conserved and the positivity of volume fractions is preserved [13].

In this chapter, we use extensive numerical examples to indicate that the CESE method captures shocks and contact discontinuities sharply without spurious oscillations and proves to be a robust and accurate numerical tool for studies of compressible multi-fluid flows.

#### **7.1 Richtmyer–Meshkov Instability**

The Richtmyer–Meshkov instability occurs when an initially perturbed interface separating two different fluids is impulsively accelerated. This impulsive acceleration, in RMI studies especially, is provided by shock waves mostly. Two basic mechanisms dominate the RMI development, are the baroclinic vorticity and the pressure disturbance. For the baroclinic vorticity, the extent of the misalignment of the pressure gradient across the shock with the density gradient across the material interface makes contribution to the perturbation growth. While for the pressure disturbance, the wave system plays an important role.

As the inertial confinement fusion (ICF) cares more about the interaction of a converging shock with a disturbed interface, the converging RMI has become an imperative [22, 23]. The nature of geometrical convergence in converging RMI, however, makes the perturbation development more complicated because of the coupling of the Bell–Plesset (BP) effect [24, 25], Rayleigh–Taylor (RT) effect [26, 27], and the multiple shock impacts therein.

In perfect agreement with their experimental images [28], Zhai et al. [17] used the upwind CESE method to make further discussion about the converging RMI issues. Figure 7.1 shows the shock behaviors and interface morphologies after a converging shock interacting with perturbed interface. The inner test gas is SF6 and the interface initial amplitude is 1 mm. The comparison clearly approved the applicability of the method used in simulating converging shock waves and converging RMI issues.

The Rayleigh–Taylor effect on the phase inversion (RTPI) was the major concern in Zhai's work [17], which is supposed to be helpful in finding a freezing state interface in ICF physics. In this work, the influence of the initial amplitude on the RTPI was firstly investigated. Figure 7.2 depicts three converging RMI cases with different initial amplitudes (Cases I1, I2 and I3 with *a*0 = 1.0, 1.65, and 2.0 mm, respectively), according to which, Fig. 7.3 records the time-variation of the interface amplitudes. It was found that the amplitude histories can be roughly segmented into three stages after the initial shock compression. At the first stage, the amplitude of the perturbed interface increases because of the RMI, while during the second stage before re-shock, the amplitude reduces owing to the RT stabilization effect which is caused by the stronger adverse pressure gradient near the geometry origin. Based

**Fig. 7.1** The evolution of a single-mode interface with initial amplitude 1 mm. upper: experimental schlieren image; lower: numerical schlieren image at the same instants. (IS = incident shock, II = initial interface, RS = reflected shock, TS = transmitted shock, SI = shocked interface, RTS = reflected transmitted shock, STS = secondary transmitted shock). Courtesy of Z. G. Zhai [17]

on the amplitude histories, they concluded that whether a normal phase inversion occurs or not depends on the competition between the RT stabilization effect and RM instability caused by the re-shock. Especially, for a special initial amplitude, there is a critical state with a zero-amplitude of the interface at the re-shock.

Furthermore, a parametric study was performed to evaluate the influences of Atwood number, shock Mach number, and initial radius of the interface on the variation of critical amplitude with the azimuthal mode number. Their thorough investigation well explained the reason why this sensitive RTPI phenomenon was not observed in previous converging RMI studies.

The even more complicated converging RMI at a dual-mode interface was numerically studied by Zhou et al. [18] using the same upwind CESE methodology. It is not surprising that the dual-mode or multi-mode interface is far more complicated than the single-mode converging counterpart due to the existence of mode interferences such as harmonic generation and bubble merger. By comparing the detailed processes of the interface deformation and wave propagation for different dual-mode cases together with the development of a referencing single-mode interface, significant influence of the phase difference between two basic waves on the instability development can be clearly distinguished. While after the re-shock, the discrepancy becomes much smaller due to the weak dependence of the re-shocked RMI on the pre-re-shock state.

**Fig. 7.2** Schlieren images showing the evolution of the single-mode interface. (PI = phase inversion, RTPI = Rayleigh–Taylor phase inversion). Courtesy of Z. G. Zhai [17]

The study on dynamics of the dual-mode converging RMI clearly reveals the mode coupling effect, in which the growth of the first mode was found to be inhibited, promoted, or not influenced, depending upon the first mode amplitude as well as the phase difference between the two basic waves. These findings, including all the other parametric studies in their work [18], were considered to be of great help for designing an optimal structure of the ICF capsule.

**Fig. 7.3** Time-variation of interface amplitude in three perturbed cases. (PI = phase inversion, RTPI = Rayleigh–Taylor phase inversion). Courtesy of Z. G. Zhai [17]

#### **7.2 Rayleigh–Taylor Instability**

Other than the shock-induced interface instability, the rarefaction-induced instability was performed by [29] using the upwind space–time CESE method. Some of their numerical schlieren images of the SF6/air interface instabilities induced by the rarefaction waves are presented in Fig. 7.4, in which the rarefaction waves were generated by simulating a diaphragm burst in a SF6 shock tube. The influence of rarefaction wave acting time was demonstrated by varying the distance between the initial diaphragm and the gas interface (shown in cases SA1, SA3, and SA5); while the influence of the rarefaction wave strength was demonstrated by setting different pressure ratios between the two sides of the initial diaphragm (shown in SA3, SA7, and SA9).

Generally, after the rarefaction wave sweeps the SF6/air interface from the downside, RTI is triggered and the perturbation on the interface gradually grows. Different from the classical RMI where an impulsive action is exerted, and also different from the classical RTI where a continuous body force exists, the rarefaction wave imposes acceleration to the interface within a limited time. When the rarefaction wave leaves the interface, the baroclinic vorticity dominates the later interface evolution and KHI occurs on the interface as the further stretching of spikes and bubbles.

Based on the numerical results provided by the upwind CESE method, theoretical models were modified in Liang's work [29] considering the time-varying acceleration and the growth rate transition from RTI to RMI. It was also found that the interface perturbation can be more unstable under the rarefaction wave condition than under

**Fig. 7.4** Schlieren images of the SF6/air interface evolution induced by rarefaction waves in cases **a** SA1, **b** SA3, **c** SA5, **d** SA7, and **e** SA9. Numbers represent the time in µs. **f** schlieren images of the interface evolution in cases SA1 (left), SA3 (middle) and SA5 (right) at dimensionless time 10.3. (RWH = rarefaction wave head, RWT = rarefaction wave tail). Courtesy of Y. Liang [29]

the shock wave condition due to the larger amount of vorticity deposited by the continuous pressure gradient.

#### **7.3 Shock Refraction**

The upwind space–time CESE method has also been used in simulating the shock refraction phenomenon at an inclined air/helium interface in the cylindrical converging shock scenario [16]. Figure 7.5 presents the experimental and numerical schlieren images of the shock refraction, in which good agreement between the two was achieved. Moreover, the numerical results apparently provide much cleaner images than their experimental counterparts (especially the region behind the deformed gas interface).

It is observed that, during the incident shock wave converging, when the shock velocity increases and the incident angle (with respect to the gas interface) decreases,

**Fig. 7.5** Sequences of experimental (lower) and numerical (upper) schlieren frames showing the evolution of a converging shock wave refracting at a 45° tilted interface. (*i* = incident shock, *m* = material interface, *t* = transmitted shock). Courtesy of Z. G. Zhai [16]

the wave pattern is found to transit from one to another. In the work of Zhai et al. [16], two interfaces with 45° and 60° initial incident angles were examined. For the 45° case, the shock pattern was found to transit from free precursor refraction (FPR) to bound precursor refraction (BPR), and then to regular refraction with reflected shock (RRR); while for the 60° case, it was found to be from twin von Neumann refraction (TNR), to twin regular refraction (TRR), to free precursor von Neumann refraction (FNR), and finally to FPR. These clearly depicted transition sequences that never occur in the planar shock scenarios, greatly enriched the shock refraction classification.

#### **7.4 Shock–Gas-Bubble Interaction**

Here, we review the numerical studies of shock-accelerated inhomogeneous flows conducted by Shen et al. [13], Fan et al. [15], and Guan et al. [14] first. Figure 7.6 shows a general schematic of the computational setting of a shock-accelerated inhomogeneity. The homogeneity can be gas bubble [13], water column [13] or droplet [14], and the shape of the initial homogeneity can be circular, square, rectangular, and even triangular [15].

**Fig. 7.6** Schematic of the initial setup for the interaction of an incident shock wave and a gas/water inhomogeneity

In the shock–helium-cylinder case [13], a planar shock propagates from left to right at Mach 1.22 and impacts the helium cylinder with a radius of 25 mm. The cylinder morphology after shock impact is shown in Fig. 7.7, in which the gas interface is sharply captured. The vortices, deposited due to baroclinicity, are clearly seen developing on the interface. As is shown, a salient feature of the helium cylinder morphology is the piercing of a middle air jet that finally separates the whole cylinder into two symmetric lobes.

A quantitative validation of the CESE method in simulating the shock–heliumcylinder interaction was performed by comparing the trajectories of the upstream point, the jet point, and the downstream point of the helium cylinder with the ones obtained by the level set method [30]. Good agreement of the results was achieved between the two methods, as depicted in Fig. 7.8.

The maximum-principle-satisfying upwind CESE scheme was applied by Fan et al. [15] to present a comprehensive study of jet-formation phenomenon in the interaction of a planar shock with a variety of heavy gas inhomogeneities. Comparing with the above light helium cylinder scenario, the heavy R22 cylinder experiences an obvious different deformation pattern. Figure 7.9 shows the cylinder morphology of the R22 cylinder with a radius of 25 mm impacted by a planar incident shock with Mach number 1.22 (to numerically reproduce the experiment conducted by Haas and Sturtevant [31]). The shocks that transmitted into this heavy cylinder were found converging at the downstream pole. The ensuing shock pattern and pressure gradient leads to a jet at the downstream cylinder surface.

Fan et al. [15] summarized that the heavy gas jet forms not only in the scenario of circular cylinder, but also in a series of differently shaped heavy gas cylinders, including square, rectangle, and triangle. Fan's square case reproduced the experiment performed by Luo et al. [32] where the shock Mach was 1.17 in air and the side length of the SF6 square was 56.6 mm. To reveal more details of the shock converging patterns in a small region, a slightly over-fined mesh was used in their simulation, which greatly reduced the numerical viscosity, making the secondary vortices at the late stage more pronounced than its corresponding experimental images. Other than that, numerical results shown in the lower halves of Fig. 7.10 agree quite well with the upper halves' experimental images, including the gas interface morphology, the shock patterns, and the consequent jet formation at the middle leeward surface of

**Fig. 7.7** Deformation history of a helium bubble impacted by a Mach 1.22 shock. Courtesy of H. Shen [13]

**Fig. 7.8** Position history of the upstream point, the jet point, and the downstream point of the helium bubble. Courtesy of H. Shen [13]

**Fig. 7.9** Shadowgraphs of the evolution of an R22 circular cylinder under the impact of an incident shock with Mach number 1.22. (*is* = incident shock; *dts* = direct transmitted shock; *rs* = reflected shock; *rts* = reversed transmitted shock; *sl* = slip line.) Courtesy of E. Fan [15]

the square. Based on the detailed information provided by the CESE simulation, the mechanism of jet formation was revealed. It was found that for all the cases they simulated, the key factor in forming a downstream jet is the formation of type II shock-shock interaction [33]. According to this, a geometrical criterion was proposed to determine whether a jet will be formed [15].

#### **7.5 Shock–Water-Droplet Interaction**

The ability of the upwind CESE method to handle gas–liquid interfaces was initially demonstrated in Shen's work [13]. The setup of this problem resembles that of the shock–bubble interaction (see Fig. 7.6). Due to the high-density ratio and large difference between the thermodynamic properties of air and water, more challenges

**Fig. 7.10** Comparison of experimental images (upper) by Luo et al. [32] and the numerical results (lower) of the planar incident shock impact on a square cylinder using CESE method. (*is* = incident shock, *rs* = refracted shock, *ts* = transmitted shock, *ms* = Mach stem, *gi* = gas interface, *dts* = direct transmitted shock, *lts* = lateral transmitted shock, *ds* = diffracted shock, *uvr* = upstream vortex pair, *dvr* = downstream vortex pair). Courtesy of E. Fan [15]

are there in simulating the shock–water-column interaction. In this case, the radius of the water column was 2.4 mm and a shock wave with Mach number 1.47 was launched to sweep over the column from left to right. Upon the interaction of a shock with water column (or water droplet), two essential issues are extensively discussed. One issue concerns the immediate interaction of the shock wave (or rarefaction wave) with the water column (or droplet). It covers an extremely short period of time compared to the whole water column breakup, in spite of its vital role in the following interface deformation. The other issue concerns the long-period water column deformation, which covers from the initial protrusion growth on the water column surface to the later large interface distortion and mass loss of the bulk of the column. At the early stage, as shown in Fig. 7.11, the shock propagates through the water column as though passing over a rigid cylinder, nearly no interface deformation is detected. However, the shock waves as well as rarefaction waves bounce back and forth inside the water column, making the internal pressure change dramatically and becomes highly heterogeneous. After the incident shock has passed, as shown in Fig. 7.12, atomized water is gradually stripped away by the high-speed post-shocked flow. The water column develops into a crescent shape and the downstream was covered up by the transversely spreading atomized water. Note that the instabilities were captured in detail in this CESE simulation. The history of drag coefficient of the water column was also derived from the numerical results, as shown in Fig. 7.13. Again, good agreement was obtained in comparison with numerical simulations of Chen [34] and experimental data of Igra and Takayama [35].

**Fig. 7.11** Sequences of pressure contours (unit: Pa) in the early stage of a shock–water-column interaction. Courtesy of H. Shen [13]

**Fig. 7.12** Evolution of the numerical schlieren plot during the stripping breakup of the water column. Courtesy of H. Shen [13]

The case of a spherical water droplet deformation under planar shock impact was studied by Guan et al. [14] using axisymmetric simulation. As shown in Fig. 7.14, the initial protrusions emerging on the droplet surface were clearly captured by the CESE methodology. The internal flow pattern of the droplet was visualized numerically in this study. It is seen in Fig. 7.15 that the internal flow pattern of the droplet is established soon after the impact by the incident shock and is held steady for a long

**Fig. 7.13** Drag coefficient *Cd* history of a water column impacted by a shock versus dimensionless time *t \**. Courtesy of H. Shen [13]

**Fig. 7.14** Comparison of the numerical (lower part) and experimental (upper part) results at different instants. (EQ = equator; WS = windward stagnation; LS = leeward stagnation; C = corrugation; KHI = Kelvin–Helmholtz instability, AT = atomization). Courtesy of B. Guan [14]

time. For the first time, a saddle point was observed in this internal flow pattern. Further, a simple theory was proposed to correlate the stationary position of the saddle point with the Mach number of the incident shock.

**Fig. 7.15** Morphologies and internal flow patterns of a water droplet under shock impact with different shock Mach numbers *Ms*. Contours show the density of air with units kg/m3. Courtesy of B. Guan [14]

#### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **Chapter 8 Application: Detonations**

Detonation is a shock-induced combustion in which chemical reactions are closely coupled with shock waves. The shock wave compresses the reactant with an abrupt increase in temperature and pressure, initiating the reactants to be burnt into products. The intense heat release permits the high propagating speed of the shock wave to be sustained. It is fundamental research related to both the safety industry and propulsion systems. For most explosive mixtures, detonation wave speeds are formulated by Chapman–Jouguet (CJ) theory. Typical detonation velocities for gaseous mixtures generally range from 1400 to 3000 m/s. Behind the shock, the time scale for reactions is commonly on the order of microseconds or even less. Furthermore, the detonation front is intrinsically unstable, forming transient multi-dimensional structures. Many studies revealed that high resolution is necessary to resolve the essential detonation structures. Due to its complex nature and multiple time scales, detonation is thus a challenging problem for solvers on shock-capturing capability, robustness, and computational efficiency. This chapter will present several essential aspects of detonation research by applying the CESE schemes.

#### **8.1 Gaseous Detonations**

In detonation simulations, the requirement for the computational resources is usually very high because the satisfactory resolution is necessary to compute the coupling between flow and chemical reactions and to resolve the detonation structures. As a result, applying detailed chemistry is usually limited to problems with simple and small domains or particular situations. Otherwise, simulations with simplified models are preferred to gain insight into the physics without loss of credibility. For illustrative purposes, we consider a system for ideal gas and a simplified chemical model. The 2D inviscid reactive compressible flow equations can be expressed as

96 8 Application: Detonations

$$\frac{\partial \mathbf{U}}{\partial t} + \frac{\partial \mathbf{F}}{\partial x} + \frac{\partial \mathbf{G}}{\partial \mathbf{y}} = \mathbf{S} \tag{8.1}$$

where the conserved variable vector *U*, the flux vectors *F*, *G* and source term vector *S* are defined as *U* = [ρ, ρ*u*, ρ*v*, *E*, ρλ] *<sup>T</sup>*, *F* = [ρ*u*, ρ*u<sup>2</sup>*+ *p*, ρ*uv*, (*E* + *p*)*u*, ρλ*u*] *T* , *G* = [ρ*v*, ρ*uv*, ρ*v*<sup>2</sup>+ *p*, (*E* + *p*)*v*, ρλ*v*] *<sup>T</sup>*, *S* = [0, 0, 0, 0, ω˙] *<sup>T</sup>*. The symbols ρ, *u*, *v*, *p*, *E*, and λ denote density, velocities in *x* and *y* directions, pressure, the total energy per unit volume, and mass fraction of the reactant, respectively. The perfect gas law was adopted here as *p* = ρ*RT*, where *R* is the gas constant and *T* is the gas temperature. The total energy is

$$E = \frac{p}{(\gamma - 1)} + \frac{1}{2}\rho(\mu^2 + \nu^2) + \rho\lambda\mathcal{Q},\tag{8.2}$$

where γ and *Q* are the specific heat ratio and the specific heat release. The chemical reaction rate is formulated by the one-step Arrhenius model as

$$
\dot{\phi} = -k\rho\lambda e^{-E\_a/RT},\tag{8.3}
$$

where *k* is the pre-exponential factor, and *Ea* denotes the activation energy.

#### *8.1.1 Ignition of Detonations*

Ignition of detonation refers to the formation of a detonation wave through either an instantaneous onset or the transition from deflagration. Shen and Parsani [1] studied the direct ignition of spherical detonation by using the upwind CESE scheme to solve Eq. (8.1). The distance between the shock wave and the position where half of the reactant is consumed can be defined as the half-reaction length. A uniform resolution with 20 meshes per half-reaction length was adopted, and the number of mesh points is 1.6 billion. It is well-known that the dynamics of detonation are sensitive to the activation energy. Thus, three activation energies representing stable (*Ea* = 15), mildly unstable (*Ea* = 27), and highly unstable (*Ea* = 50) detonations were studied. For different *Ea*, *k* was adjusted to fix the half-reaction length in unit length. A hot spot with *ps* = (γ − 1)*Es*, *Ts* = 20.0, and a radius of 1 was used to form a blast wave to ignite the mixture. All variables above are nondimensionalized with respect to the state of the unburnt reactant.

The peak pressure history behind the shock for 1D detonations are depicted in Fig. 8.1, and typical 2D contours are plotted in Fig. 8.2. For stable case (Fig. 8.1a), the subcritical, critical, and supercritical regimes are observed successively when *Es* increases from 1.8 × 104 to 1.0 × 105. In the critical regime, *Es* = 2.0 × 104 (Fig. 8.1b), *P*sh first decays and then rapidly increases due to the formation and amplification of the pressure pulse behind the leading shock. The 2D result shows a

**Fig. 8.1** Peak pressure history at the shock front (*P*sh) as a function of position and source energy (*Es*), for **a** *Ea* <sup>=</sup> 15, **b** *Ea* <sup>=</sup> 15, *Es* <sup>=</sup> <sup>2</sup><sup>×</sup> 104, **c** *Ea* <sup>=</sup> 27, **d** *Ea* <sup>=</sup> 50. Courtesy of H. Shen [1]

highly oscillatory pattern for *P*sh due to the development of multidimensional instabilities, and the run-up distance deviates from 1D solution. Three critical energies were detected for the mildly unstable case (Fig. 8.1c), The authors further observed the fourth critical energy when refining the search using a minor incremental *E*s. However, there might be a unique critical energy in the 2D case. The inconsistency between 1 and 2D becomes more apparent for the highly unstable cases. For all the *Es*  tested, 1D detonation eventually failed (Fig. 8.1d). The results suggest that, without large overdriven, the 1D detonation propagation through auto-ignition is impossible for the highly unstable cases. Contrarily, the critical energy is approximately 1.5 × 105 for 2D cases (Fig. 8.2d). The key factors dominating the direct ignition of 2D detonation are the unsteadiness arising from the decay of the leading shock, heat release from the chemical reaction, and the inherent multidimensional instabilities. In 1D detonations, only the first two aforementioned factors exist. Detonation fails when the excessive unsteadiness overtakes the heat release. The transverse waves induced by the multidimensional instabilities create local over-driven detonations. On the one hand, these local over-driven detonations facilitate the propagation of the global detonation; on the other hand, induce stronger expansion waves to increase the risk of failure. The competition among these three factors should be considered to provide a comprehensive model of direct initiations.

**Fig. 8.2** Maximum pressure history in two dimensions. **a** *Ea* <sup>=</sup> 15, *Es* <sup>=</sup> 2.5 <sup>×</sup> 104, **b** *Ea* <sup>=</sup> 27, *Es* <sup>=</sup> 4.5 <sup>×</sup> 104, **c** *Ea* <sup>=</sup> 50, *Es* <sup>=</sup> 1.0 <sup>×</sup> 104, **d** *Ea* <sup>=</sup> 50, *Es* <sup>=</sup> 1.5 <sup>×</sup> 104. Courtesy of H. Shen [1]

Another situation of detonation ignition occurs when the gaseous reactant in the unconfined space is ignited through the detonation wave emerging from a small tube. When the detonation wave diffracts from a tube into the unconfined space, the substantial expansion will perturb the detonation front. If tube size is small enough, the detonation wave will quench (subcritical), otherwise, the detonation will continue to propagate (supercritical). Shi et al. [2] investigated the 2D detonation diffraction using the *a*–α CESE scheme. Two scenarios were considered. If the detonation wave inside the tube does not contain transverse waves, it is classified as the planar case. Otherwise, it is classified as the cellular case. In these simulations, Eqs. (8.1)–(8.3) were employed with *Ea* = 24 to assure that 1D detonation wave is stable (For *Ea*  higher than 25.27, 1D detonation is intrinsically unstable). A resolution with 24 meshes per half-reaction length was tested converged for current settings, and the total mesh points are around 410 million. A first glance at the results revealed that the critical channel width for cellular scenarios is smaller than that for the planar case (Table 8.1). Four cases, including supercritical and subcritical for planar and cellular

#### 8.1 Gaseous Detonations 99


scenarios, were carefully discussed by examining wave structures and Lagrangian particles that are initially scattered in the flow.

For planar scenarios, the disturbance originates from the corner propagates towards the symmetric line. The diffracted shock wave decelerates and decouples with the flame front. In the shocked reactant, a compression wave is formed and amplified through the temperature gradient towards the leading shock. The compression wave is indicated by the successive abrupt change in temperature for particles E7–E9 (Fig. 8.3). When this compression wave merges with the leading shock, detonation may be facilitated if the strengthened shock wave exceeds a critical value. On the other hand, the aforementioned disturbance reflects from the symmetric line, and the reflected rarefaction (Fig. 8.4) wave may simultaneously affect the shocked reactant and hinder the formation and amplification of the compression wave. For the sub-critical case, i.e., *w* = 100, re-initiation is reproduced if the boundary condition is modified such that the rarefaction wave does not reflect towards the leading shock. Therefore, for the planar cases, re-initiation is attributed to the competition between the coalescence of the amplified compression wave with the leading shock and the strength of the reflected rarefaction wave.

For cellular scenarios, the re-initiation patterns are different from the planar cases. Despite the number of transverse waves decreases when the detonation wave expands into the unconfined space. It is interesting to note that from Fig. 8.5, a complete decoupling is not observed during the propagation of the supercritical case (*w* = 85) (a complete decouple is observed before re-initiation in the planar *w* = 110

**Fig. 8.3** Temperature traces of particles along the re-initiation direction for planar detonation diffraction: **a** *w* = 110 and **b** *w* = 100. Courtesy of L. S. Shi [2]

**Fig. 8.4** Numerical schlieren images for planar case with *w* = 100. Red dash lines indicated the head of rarefaction. Courtesy of L. S. Shi [2]

**Fig. 8.5** Numerical soot foils for **a** planar detonation diffractions, *w* = 110, **b** cellular detonation diffractions, *w* = 85. The red dashed lines indicate the disturbance line predicted by Skews' model. The blue dash-dot line is the re-initiation path along which particles *Ei* are initially located. Courtesy of L. S. Shi [2]

case). Detonation wave sustains through the continuous formation of local overdriven detonation when transverse waves collide. Previous experimental research revealed that the critical tube diameter for unstable mixtures is *dc* ≈ 13λ (λ: detonation cell size), but *dc* ≈ 30λ or even more were measured for stable mixtures highly diluted with argon. The distinct difference between the re-initiation mechanisms planar and cellular detonation diffractions elaborated in current simulation provides a plausible explanation of the discrepancy of the correlations on critical tube sizes. In unstable detonations, where the formation of local hot spots is essential to sustain the detonation, the correlation between the *dc* and λ is established. It is well-known that the cellular patterns in stable detonation are regular, with transverse wave velocity approximately being the sound velocity in the burned products. This corresponds to the planar scenarios. Detonation re-initiates if the influence of the reflected rarefaction wave is negligible. Thus the *dc* / λ correlation may not be necessary.

#### *8.1.2 Dynamics of Detonations*

As discussed in Sect. 8.1, multi-dimensional instability is a fundamental phenomenon in detonation. In this section, we will elaborate more on the studies on detonation propagation in one, two, and three dimensions using CESE schemes.

Wang et al. [3] have implemented 2D CESE on detonation reflection on a wedge using a detailed reaction model. Later, in the study on the propagation modes of stoichiometric H2/O2 3D cellular detonation in a square duct [4], two types of propagation modes were verified, which depended on the configuration of the initial conditions. Unreacted pockets with complex structures were also captured by using the improved CESE scheme.

Numerous numerical studies based on detailed chemical mechanisms have quantitatively reproduced the key structures of multi-dimensional detonations. They experienced barriers to predict the cell size correctly. For example, for H2/O2 mixture highly diluted by argon, the regular cell sizes predicted from simulation using detailed chemical mechanisms were approximately half of those measured in experiments. One of the important physical phenomena in high-speed flow, the vibrational nonequilibrium of molecules, was usually ignored in early numerical studies. Behind the detonation wave, the time scales for vibrational relaxation and chemical reaction are comparable. Shi et al. [5] re-examined this problem with the consideration of vibrational relaxation and its coupling with reactions. The convection of vibrational energy was integrated into governing equations, and the translational-rotational energy and the vibrational energy were separated in the formulation of the total energy, which is composed of enthalpy of formation, translational-rotational energy, kinetic energy, and vibrational energy. The energy exchange rate was calculated by the Landau-Teller model.

Four scenarios were compared in 1D and 2D simulations. (1) The mixture is assumed to be thermodynamically equilibrium. (2) Vibrational relaxation is incorporated, and the translational temperature is kept as the dominant temperature of the chemical reactions. (3) Vibrational relaxation is considered, and the geometric average temperature is used to compute the reaction rates. (4) The physically consistent vibration-chemistry-vibration coupling model is adopted to account for the effect of vibrational non-equilibrium on chemical reaction rates. Examples of 1D detonation structures are depicted in Fig. 8.6. Model (2) exhibits a temperature overshot behind the shock compared to the equilibrium case because the mixture remains vibrational cold. The vibrational temperature approaches the equilibrium state gradually, and a noticeable disparity still exists when exothermic reactions are triggered. Since translational-rotational temperature controls reactions in model (2), the overshot in temperature results in a marginal decrease in the half-reaction length (δ). In model (3), the geometric averaged temperature is used to calculate the reaction rate. The onset of severe chemical reactions occurs only when the vibrational relaxation approaches equilibrium. The half-reaction length increases by 1.64 times that in the model (1). A similar trend is observed in the case using model (4). The detonation

**Fig. 8.6** Temperature and H2 distribution with initial condition of 0.1 atm and 300 K. Left: scenarios 1 and 2. Right: scenarios 1 and 3. Courtesy of L. S. Shi [5]

cell widths are predicted from 2D simulation using models (1) and (2) are approximately the same. The cell widths using models (3) and (4) are approximately 1.32 and 1.34 times the value of that in the model (1). It reveals that the vibrational-chemical coupling effect may be one of the factors responsible for the failure to correctly predict the cell size.

Later, Uy et al. [6] further examine the effect of vibrational non-equilibrium on the one-dimensional instabilities using the CESE scheme. In this study, 1D pistonsupported detonation case with fixed non-dimensional specific heat release *Q* = 50, the ratio of specific heats γ = 1.2, and characteristic vibrational temperature θ = 20. The ratio between the chemical time scale τ *<sup>c</sup>*and the vibrational time scale τ *<sup>V</sup>*was selected as a control parameter. Park's two-temperature model was used to evaluate the effect of vibrational non-equilibrium on the chemical reaction rate. The neutral stability limit under the vibrational equilibrium assumption is approximately *Ea* = 26.47, above which the detonation is unstable with amplified oscillation. Numerical tests reveal *Ea* = 27 with the thermal equilibrium assumption being longitudinally unstable. Nevertheless, for finite vibrational relaxation rate (i.e., the τ <sup>α</sup> ≡ τ*c*/τ<sup>V</sup> /= ∞), the amplitudes of pressure oscillations decay for smaller τ <sup>α</sup> = 3, 5, 7, and the oscillation decays much faster for as τ <sup>α</sup> decreases. The results imply that the propagation of detonation is stabilized because of the vibrational nonequilibrium. Theoretical analysis using linear stability analysis (LSA) under different overdriven factors*f* =(*D*/*D*CJ) 2 was also provided in Uy et al. [7]. The neutral stability limit predicted using the CESE scheme agrees very well with LSA (Table 8.2). The accuracy of CESE in solving detonation problems is further confirmed.

**Table 8.2** Comparison of the neutral stability limit (NSL) and the period of oscillation (PO) computed by LSA and numerical simulation. Case I: *f* = 1, thermal equilibrium (eq). Case II: *<sup>f</sup>*<sup>=</sup> 1, thermal nonequilibrium (neq), <sup>τ</sup> <sup>α</sup> <sup>=</sup> 5. Case III: *f* <sup>=</sup> 1, neq, <sup>τ</sup> <sup>α</sup> <sup>=</sup> 7. Case IV: *f* <sup>=</sup> 1, neq, <sup>τ</sup> <sup>α</sup> <sup>=</sup> 9. Case V: *Ea* <sup>=</sup> 50, eq. Case VI: *Ea* <sup>=</sup> 50, neq, <sup>τ</sup> <sup>α</sup> <sup>=</sup> 5. Case VII: *Ea* <sup>=</sup> 50, neq, <sup>τ</sup> <sup>α</sup> <sup>=</sup> 10. Case VIII: Ea <sup>=</sup> 50, neq, <sup>τ</sup> <sup>α</sup> <sup>=</sup> 20. Courtesy of C. K. Uy [7]


#### *8.1.3 Rotating Detonation Waves*

Because of the rapid compression of the mixture and the quasi-constant volume combustion, the high thermodynamic efficiency of detonation makes it favourable to be employed in developing potential detonation engines. The rotating detonation engine (RDE) is one of the promising candidates in which the detonation wave propagates circumferentially, and the reactant is continuously injected into the combustion chamber.

Figure 8.7 depicts the flow features of RDE simulated with the upwind CESE schemes. The 2D simulation is essentially the unwrapped RDE channel from a 3D geometry. The top is the RDE outlet described by a non-reflective boundary. The reactant is injected into the combustor from the bottom through micro-Laval nozzles. The ideal injection was assumed that all the bottom areas could inject fresh reactant. Typical flow features of RDE can be observed, characterized by a rotating detonation wave, triangular reactant layer, oblique shock, and slip line. Furthermore, Wang et al. [8] used an improved CNI 2D CESE scheme to study the kerosene/air RDE. The results suggested that decreasing the inlet/wall ratio would reduce the strength of the detonation wave, and the reaction zone would be elongated. With further decreasing the air ratio, the burned gas emerges in the triangular zone, and the detonation wave fails to be self-sustained.

**Fig. 8.7** Density contours of 2D and 3D RDE simulations using the upwind CESE scheme

#### **8.2 Two-Phase Detonations**

Liquid fuels and energetic powders are of advantages of high energy density and easy storage, which makes them common in practical propulsion systems. Detonationbased propulsion systems using these non-gaseous fuels involve reactive multi-phase high-speed flows. Although gaseous, liquid and/or solid phases are involved in these detonation phenomena, it is always referred to as two-phase detonation from the viewpoint of reactants, namely gas–liquid two-phase detonation and gas–solid two-phase detonation. Compared to gaseous detonations, two-phase detonations are characterized not only by shock waves and fast combustion, but also by multi-phase interaction and multi-scales, leading to great difficulties in numerical simulations. With the successful applications of the CESE method in gaseous detonation simulations, it is imperative to extend the CESE method to two-phase detonation problems and test its capabilities of accuracy and robustness in two-phase detonation simulations.

There are mainly two kinds of frameworks that are used to address the discrete phase in two-phase detonation simulations, namely the Eulerian–Eulerian framework and the Eulerian–Lagrangian framework. Between these two frameworks, the Eulerian–Eulerian framework, which is also referred to as the two-fluid model/method, is more common and easier in implementation and extension in two-phase detonation simulations. The discrete phase (the liquid phase in gas–liquid two-phase detonations or the solid phase in gas–solid two-phase detonations) is considered as a special continuum, such that continuum mechanics can be employed to describe the bulk motion of the discrete phase. Then, the discrete phase becomes another "fluid", and the gaseous-phase flow and the discrete-phase flow can be solved by similar approaches. As for the Eulerian–Lagrangian framework, the gaseous phase is solved as in the Eulerian–Eulerian framework, but every discrete liquid particle (droplet) or solid particle is tracked individually by Newton's laws of motion.

In two-phase detonation modelling, it is a general way to assume that the particle (liquid or solid) of the discrete phase is small, spherical in shape and uniformly suspended in the gas atmosphere. Additionally, the two-phase suspension is assumed to be diluted enough to neglect the volume fraction of discrete particles and particle–particle collisions. Two-phase detonation is a highly transient two-phase flow problem involving strong shock waves; hence, thermal and mechanical non-equilibrium between the gas and particles should be considered. Further, uniform distribution of temperature is considered within particles due to their small particle sizes, and ideally, all reaction heat is absorbed by gas only. Based on the above assumptions, the two-dimensional governing equation of the gaseous phase in an Eulerian–Eulerian framework can be expressed as follows [9]:

$$\frac{\partial \mathbf{U}}{\partial t} + \frac{\partial \mathbf{F}}{\partial \mathbf{x}} + \frac{\partial \mathbf{G}}{\partial \mathbf{y}} = \mathbf{S} + \mathbf{W},\tag{8.4}$$

where *U* is the vector of conserved variables, *F* and *G* the conservation flux vectors in the *x*- and *y*-directions, *S* the vector of two-phase interaction source terms, and *W*  the vector of chemical reaction source terms, respectively. They can be expressed as

$$\mathbf{U} = \begin{bmatrix} \rho\_1 \\ \dots \\ \rho\_{ns} \\ \rho u \\ \rho v \\ E \end{bmatrix}, \quad \mathbf{F} = \begin{bmatrix} \rho\_1 u \\ \dots \\ \rho\_{ns} u \\ \rho u^2 + p \\ \rho uv \\ (E+p)u \end{bmatrix}, \quad \mathbf{G} = \begin{bmatrix} \rho\_1 v \\ \dots \\ \rho\_{ns} v \\ \rho uv \\ \rho v^2 + p \\ (E+p)v \end{bmatrix}, \quad \mathbf{W} = \begin{bmatrix} \dot{\boldsymbol{\alpha}}\_1 \\ \dots \\ \dot{\boldsymbol{\alpha}}\_{ns} \\ 0 \\ 0 \\ 0 \end{bmatrix}, \tag{8.5}$$

and

$$S = \begin{bmatrix} 0 \\ \vdots \\ J\_p \\ \vdots \\ 0 \\ -f\_x + u\_p J\_p \\ -f\_y + v\_p J\_p \\ -q\_p - (u\_p f\_x + v\_p f\_y) + (E\_p/\rho\_p) J\_p \end{bmatrix}. \tag{8.6}$$

In the above equations, ρ*i* is the mass density of species *i* and *i* = 1, …, *ns*; *ns* is the number of species contained in the gas mixture; *p*, *u*, and *v* are the gas pressure and *x*- and *y*-components of gas velocity, respectively; ω˙*i* is the mass production rate of gaseous species *i* by chemical reactions; and ρ and *E* are the mass density and total energy per unit volume of the gas mixture, expressed as

106 8 Application: Detonations

$$\rho = \sum\_{i=1}^{ns} \rho\_i, \quad E = \rho h - p + \frac{1}{2}\rho \left(\mu^2 + \nu^2\right), \tag{8.7}$$

where *h* is the specific enthalpy of the gas mixture, calculated by

$$h = \sum\_{i=1}^{ns} \frac{\rho\_i}{\rho} h\_i,\tag{8.8}$$

with the specific enthalpy of each individual species *hi* as a function of gas temperature *T*, obtained from a species thermodynamic data base, such as the NASA Glenn data base [10]. By further assuming each species performs as a perfect gas, the equation of state of the gas mixture is then given by

$$p = \sum\_{i=1}^{ns} \rho\_i R\_i T, \quad R\_i = \frac{R\_0}{W\_i} \tag{8.9}$$

where *R*0 = 8.314 J/(mol K) is the universal gas constant, and *Wi* is the molar mass of species *i*. Variables in Eq. (8.6) are related to the properties of the discrete phase and the interaction between two phases. Their definitions will be given in the following paragraphs.

As for the discrete phase, the governing equation has a similar form as Eq. (8.4):

$$\frac{\partial \mathbf{U}\_p}{\partial t} + \frac{\partial \mathbf{F}\_p}{\partial \mathbf{x}} + \frac{\partial \mathbf{G}\_p}{\partial \mathbf{y}} = \mathbf{S}\_p,\tag{8.10}$$

where *Up*, *Fp* and *Gp* are the vectors of conserved variables, conservation fluxes in the *x*- and *y*-directions of the discrete phase, which can be given by

$$\begin{aligned} \mathbf{U}\_{P} &= \begin{bmatrix} \rho\_{p} \\ \rho\_{p}u\_{p} \\ \rho\_{p}v\_{p} \\ E\_{p} \\ N\_{p} \end{bmatrix}, \quad \mathbf{F}\_{P} = \begin{bmatrix} \rho\_{p}u\_{p} \\ \rho\_{p}u\_{p}^{2} \\ \rho\_{p}u\_{p}v\_{p} \\ E\_{p}u\_{p} \\ N\_{p}u\_{p} \end{bmatrix}, \quad \mathbf{G}\_{P} = \begin{bmatrix} \rho\_{p}v\_{p} \\ \rho\_{p}u\_{p}v\_{p} \\ \rho\_{p}v\_{p}^{2} \\ E\_{p}v\_{p} \\ N\_{p}v\_{p} \end{bmatrix}, \\ \mathbf{S}\_{P} &= \begin{bmatrix} -J\_{p} \\ f\_{x} - u\_{p}J\_{p} \\ f\_{y} - v\_{p}J\_{p} \\ q\_{p} + \left(u\_{p}f\_{x} + v\_{p}f\_{y}\right) - \left(E\_{p}/\rho\_{p}\right)J\_{p} \\ 0 \end{bmatrix} \end{aligned} \tag{8.11}$$

In Eqs. (8.6) and (8.11), *Jp* is the mass regression rate (the combustion rate) of the discrete phase and it is determined by the combustion model of the discrete phase. ρ*p* is the density of the discrete phase in the suspension. *up* and *vp* are the *x*- and

*y*-components of the velocity of the discrete phase, respectively. *Ep* and *Np* are the total energy per unit volume and the particle number density of the discrete phase, respectively. They can be calculated by

$$E\_p = \rho\_p e\_p + \frac{1}{2}\rho\_p (\mu\_p^2 + \nu\_p^2), \quad N\_p = \frac{6\rho\_p}{\pi \rho\_m d\_p^3},\tag{8.12}$$

where *ep* is the specific internal energy of the discrete phase and can also be obtained from a species thermodynamic data base as a function of discrete phase temperature *Tp*; *dp* is the diameter of the particle; and ρ*m* is the material density of the discrete phase.

Additionally, the *x*- and *y*-components of drag force acting on the discrete phase, *f x* and *f y*, can be modelled as follows:

$$\begin{cases} f\_x = N\_p \frac{\pi}{8} C\_D d\_p^2 \rho \left| V - V\_p \right| \left( \mu - u\_p \right) \\\ f\_y = N\_p \frac{\pi}{8} C\_D d\_p^2 \rho \left| V - V\_p \right| \left( \nu - v\_p \right) \end{cases} \tag{8.13}$$

where *CD* is the drag coefficient,

$$C\_D = \begin{cases} \frac{24}{\text{Re}\_p} \left( 1 + \frac{1}{6} \text{Re}\_p^{2/3} \right), \text{ for } \text{Re}\_p < 1000\\ 0.424, & \text{for } \text{Re}\_p \ge 1000 \end{cases} \tag{8.14}$$

In Eqs. (8.13) and (8.14), the relative velocity and the relative Reynolds number Re*p* between the gas and discrete phases can be calculated by

$$\left|\boldsymbol{V} - \boldsymbol{V}\_{p}\right| = \left[\left(\boldsymbol{u} - \boldsymbol{u}\_{p}\right)^{2} + \left(\boldsymbol{v} - \boldsymbol{v}\_{p}\right)^{2}\right]^{1/2},\tag{8.15}$$

and

$$\text{Re}\_p = \frac{\rho d\_p \left| V - V\_p \right|}{\mu},\tag{8.16}$$

where μ is the viscosity coefficient of gas. Accordingly, the convection heat transfer between the two phases is expressed as follows:

$$q\_p = N\_p \pi d\_p \lambda \text{Nu}\_p \left( T - T\_p \right), \tag{8.17}$$

with the two-phase Nusselt number expressed as functions of Re*p* and Prandtl number Pr:

$$\text{Nu}\_p = 2 + 0.459 \,\text{Re}\_p^{0.55} \,\text{Pr} \,\text{.}\tag{8.18}$$

The CESE method has been proven to be of high accuracy and good stability in solving gaseous detonation problems in Sect. 8.1, and it is chosen to solve the twophase detonation problems as well. That is, Eqs. (8.4) and (8.10) in the Eulerian– Eulerian framework can be solved using a CESE method, such as the *a*-α scheme. Theoretically, the source terms in the governing equations can be addressed together with the space–time integration in the CESE method [11–15]. However, a separated treatment of source terms is always employed in solving two-phase detonations, because the Jacobian matrixes of the interphase interaction and chemical reaction source terms in coupling treatment are rather complicated to calculate. Notably, with the source terms treated separately, it has been proven that the good accuracy and stability of the CESE method are preserved in high-speed reactive flow simulations [5, 16, 17]. On the other hand, the characteristic time scales of the interphase interaction and chemical reaction source terms are much smaller than that of flow dynamics, always leading to stiffness problems in two-phase detonation simulations. To overcome this problem, the operator-splitting technique with multiple sub-time steps [18] is always employed and then the source terms of interphase interactions and chemical reactions are explicitly integrated as ordinary differential equations. The detailed implementation process under the Eulerian–Eulerian framework can be illustrated as follows:

$$\begin{aligned} \begin{cases} \mathbf{U}\_{\text{n}} \xrightarrow[\text{S-W}]{\text{CES}} \tilde{\mathbf{U}}\_{\text{n}+1} \\ \mathbf{U}\_{p} \xrightarrow[\text{S}\_{p}=0]{\text{CES}} \tilde{\mathbf{U}}\_{p,n+1} \end{cases} \Rightarrow \begin{cases} \Delta t' = \Delta t / \mathbf{N} \\ \mathbf{U}\_{n+1}^{(0)} = \tilde{\mathbf{U}}\_{n+1} \\ \mathbf{U}\_{p,n+1}^{(0)} = \tilde{\mathbf{U}}\_{p,n+1} \end{cases} \\\Rightarrow \begin{cases} \mathbf{U}\_{n+1}^{(\text{m})}, \mathbf{U}\_{p,n+1}^{(\text{m})} \to \mathbf{S}^{(\text{m})}, \mathbf{W}^{(\text{m})}, \mathbf{S}\_{p}^{(\text{m})} \\ \mathbf{U}\_{n+1}^{(\text{m})} = \mathbf{U}\_{n+1}^{(\text{m})} + \Delta t' [\mathbf{S}^{(\text{m})} + \mathbf{W}^{(\text{m})}] \\ \mathbf{U}\_{p,n+1}^{(\text{m})} = \mathbf{U}\_{p,n+1}^{(\text{m})} + \Delta t' \cdot \mathbf{S}\_{p}^{(\text{m})} \end{cases} \\\Rightarrow \begin{cases} \mathbf{U}\_{n+1} = \mathbf{U}\_{n+1}^{(\text{N})} \\ \mathbf{U}\_{p,n+1} = \mathbf{U}\_{p,n+1}^{(\text{N})} \end{cases}, \tag{8.19}$$

where the subscripts n and m refer to the global time step (Δ*t*) and the sub-time step (Δ*t* ' ), respectively, and N is the total number of sub-time steps within one global convection time step of the CESE method. Depending on the degree of stiffness in the problem, N can be chosen to be 10−20.

To test the accuracy and robustness of the CESE method in two-phase detonation simulations under the Eulerian–Eulerian framework, Wang et al. [15] applied the CESE method to simulate the detonation synthesis of titania (TiO2) nanoparticles. The corresponding chemical reaction can be described as follows:

$$\text{TiCl}\_4(\text{g}) + \text{O}\_2(\text{g}) + 2\text{H}\_2(\text{g}) \rightarrow \text{TiO}\_2(\text{s}) + 4\text{HCl}(\text{g})\tag{8.20}$$

As seen, it involves TiO2 solid particles as well as TiCl4, O2, H2 and HCl gases, implying that it is a gas–solid two-phase detonation problem. Simulation results showed that the simulated detonation profiles of gas density, pressure and temperature

**Fig. 8.8** Profiles of gas density, velocity, pressure and temperature in the detonation synthesis of nanosized TiO2 particles. Courtesy of G. Wang [15]

are similar to those in the gaseous detonation wave, but a sudden drop of gas density can be captured after the detonation front because of phase transition from gas to solid particles, as shown in Fig. 8.8. Additionally, it is also found that the simulated particle size of the produced TiO2 and the simulated peak pressure agreed well with other calculation and experiment data [19], respectively.

As for gas–liquid two-phase detonation, Wang et al. [3] solved the two-phase planar detonations of liquid hydrocarbon fuels with the same CESE method under the Eulerian–Eulerian framework. It is also found that the detonation profiles are similar to those of gaseous detonations, but with higher gas density, pressure and temperature behind the detonation front due to higher energy density of liquid fuels. Additionally, the length of the reaction zone is found to be longer as well, which is relative to the slower combustion rate of the liquid particle. Table 8.3 shows the comparisons of simulated detonation speeds to the experimental values. It can be revealed that the relative errors are less than 10% for particle sizes small than about 145 µm. As for two cases with larger particle sizes, the cases with relative errors exceeding 10%, it may be caused by the neglect of deformation of large liquid particles in modelling, since the combustion rate of liquid particle is influenced by particle deformation. All the above numerical results showed that that the CESE method is of high accuracy and good robustness in two-phase detonation simulations under an Eulerian–Eulerian framework.

The Eulerian–Eulerian framework is a simple and effective way to deal with the discrete phase in two-phase detonations, but it has some inherent limitations in modelling realistic two-phase suspension in industries or experiments with a specific particle size distribution, where a relatively wide range of particle diameters are involved [20–22]. Under the Eulerian–Eulerian framework, all particles within one numerical mesh are assumed to be in the same states, such as the same particle size, temperature, velocity and so on. However, the number of particles within one mesh may be large, and their states may differ depending on their initial states and interaction histories with gas phase. Moreover, the forces and the heat transfers between the gas and particles also differ by the particle size, resulting in different temperatures and


**Table 8.3** Comparisons of detonation speeds of liquid hydrocarbon fuels between simulations and experiments. Courtesy of G. Wang [3]

velocities of particles within one computational mesh. Further, it had been demonstrated that many two-phase detonation characteristics are significantly influenced by the particle size. Consequently, the Eulerian–Eulerian framework is insufficient to reflect the true physics of realistic two-phase detonations with particle size distributions and to simulate them accurately, which raises the demand of modelling realistic two-phase detonation under the Eulerian–Lagrangian framework where the discrete particles are tracked individually.

Under the Eulerian–Lagrangian framework, Eq. (8.4) is still applied as the governing equation of the gas phase, and **U**, **F**, **G** and **W** have the same forms as Eq. (8.5). However, the two-phase interaction source term **S** is expressed in the following form instead,

$$\mathbf{S} = \frac{1}{dV} \sum\_{k=1}^{np} \begin{bmatrix} 0 \\ \vdots \\ J\_{pk} \\ \vdots \\ 0 \\ 0 \\ -f\_{xk} + u\_{pk}J\_{pk} \\ -f\_{yk} + v\_{pk}J\_{pk} \\ -q\_{pk} - \left(u\_{pk}f\_{xk} + v\_{pk}f\_{yk}\right) + \frac{1}{2}\left(u\_{pk}^2 + v\_{pk}^2\right) \cdot J\_{pk} + e\_{pk} \cdot J\_{pk} \\ \left(8.214\right) \end{bmatrix},\tag{8.21}$$

where subscript *k* represents all the quantities related to the *k*th particle (solid or liquid). To include all effects of particles into the source term **S** of the gaseous equation, the summation is done within the gaseous mesh element *dV*, and *np* is the number of particles in *dV*.

The motion of every particle is then descripted using Newton's laws of motion instead of the Eulerian form Eq. (8.10). For the *k*th particle, the corresponding governing equation can be written as

#### 8.2 Two-Phase Detonations 111

$$\frac{d\mathbf{L}\_{pk}}{dt} = \mathbf{S}\_{pk},\tag{8.22}$$

where **L***pk* and **S***pk* are the vectors of the Lagrangian variables of the *k*th particle and the corresponding source terms, respectively, and they are expressed as

$$\mathbf{L}\_{pk} = \begin{bmatrix} m\_{pk}, \mathbf{x}\_{pk}, \mathbf{y}\_{pk}, m\_{pk}\boldsymbol{\mu}\_{pk}, m\_{pk}\boldsymbol{\nu}\_{pk}, \boldsymbol{E}\_{pk} \end{bmatrix}^{\mathrm{T}},\tag{8.23}$$

and

$$\mathbf{S}\_{pk} = \begin{bmatrix} -J\_{pk}, u\_{pk}, v\_{pk}, f\_{xk}, f\_{yk}, -e\_{pk}J\_{pk} + q\_{pk} \end{bmatrix}^{\mathrm{T}}.\tag{8.24}$$

In Eq. (8.23), *mpk* and *Epk* are the mass and total internal energy of the *k*th particle, respectively, and *Epk* = *mpk* ·*epk* .

Further, the Lagrangian forms of the drag force and the convection heat of the *k*th particle can be easily derived from the according Eulerian forms Eqs. (8.13) and (8.17), as follows:

$$\begin{cases} f\_{\mathbf{x}k} = \frac{\pi}{8} C\_{Dk} d\_{pk}^2 \rho \left| V - V\_{pk} \right| \left( \mu - u\_{pk} \right) \\\ f\_{\mathbf{y}k} = \frac{\pi}{8} C\_{Dk} d\_{pk}^2 \rho \left| V - V\_{pk} \right| \left( \mathbf{v} - \mathbf{v}\_{pk} \right) \end{cases},\tag{8.25}$$

and

$$q\_{pk} = \pi d\_{pk} \lambda \text{Nu}\_{pk} \{T - T\_{pk}\}. \tag{8.26}$$

Under the Eulerian–Lagrangian framework, the CESE method is again applied to solve the gas phase equation, while the source terms of interphase interactions and chemical reactions are treated separately and integrated explicitly as ordinary differential equations, along with the integration of particle Lagrangian equations, by using the operator-splitting technique as well, as depicted below,

$$\begin{aligned} \mathbf{U}\_{\mathbf{n}} \xrightarrow[\mathbf{S}=\mathbf{W}=\mathbf{0}]{\mathbf{C}\mathbf{S}\mathbf{E}} \mathbf{\widetilde{U}}\_{\mathbf{n}+1} &\implies \begin{cases} \Delta t' = \Delta t/\mathbf{N} \\ \mathbf{U}\_{\mathbf{n}+1}^{(0)} = \mathbf{\widetilde{U}}\_{\mathbf{n}+1} \\ \mathbf{L}\_{pk,\mathbf{n}+1}^{(0)} = \mathbf{L}\_{pk,\mathbf{n}} \end{cases} \\ &\Rightarrow \begin{cases} \mathbf{U}\_{\mathbf{n}+1}^{(\mathbf{m})}, \mathbf{L}\_{pk,\mathbf{n}+1}^{(\mathbf{m})} \to \mathbf{S}^{(\mathbf{m})}, \mathbf{W}^{(\mathbf{m})}, \mathbf{S}\_{pk}^{(\mathbf{m})} \\ \mathbf{U}\_{\mathbf{n}+1}^{(\mathbf{m}+1)} = \mathbf{U}\_{\mathbf{n}+1}^{(\mathbf{m})} + \Delta t' [\mathbf{S}^{(\mathbf{m})} + \mathbf{W}^{(\mathbf{m})}] \\ \mathbf{L}\_{pk,\mathbf{n}+1}^{(\mathbf{m}+1)} = \mathbf{L}\_{pk,\mathbf{n}+1}^{(\mathbf{m})} + \Delta t' \cdot \mathbf{S}\_{pk}^{(\mathbf{m})} \end{cases} \\ &\Rightarrow \begin{cases} \mathbf{U}\_{\mathbf{n}+1} = \mathbf{U}\_{\mathbf{n}+1}^{(\mathbf{N})} \\ \mathbf{L}\_{pk,\mathbf{n}+1} = \mathbf{L}\_{pk,\mathbf{n}+1}^{(\mathbf{N})} \end{cases} . \end{aligned} \tag{8.27}$$

Notably, the Eulerian–Lagrangian framework had rarely been developed to solve high-speed reactive two-phase flows, mainly because the number of fine particles has proven too large and simulations too expensive to achieve in the past. Nowadays, with the rapid development of computer technologies, parallel computation techniques such as the Message Passing Interface (MPI) technique may help to make simulations of two-phase detonations with fine particles under an Eulerian– Lagrangian framework possible. However, the application of MPI parallel computation technique to an Eulerian–Lagrangian framework is not as straightforward as that under a purely Eulerian framework. A large amount of information exchange from one computational core to another is needed for calculations of the interaction source terms between two phases, leading to formidable communication cost even exceeding the computational cost, especially when the Eulerian framework is staggered with the Lagrangian framework due to the relative movement between these two phases.

To solve the communication problem of parallel computing technique of highspeed two-phase flows under an Eulerian–Lagrangian framework, the traditional static data structure, such as multi-dimensional array, should be avoided using to store the information of the discrete phase described under the Lagrangian framework. Notably, the order of particles presented under the Lagrangian framework is not important and does not need to be preserved as its initial order. The only operation required for coding is to traverse every particle one by one. Inspired by the above facts, dynamic data structures can be introduced to store particle information and solve the communication problem under the Eulerian–Lagrangian framework. For example, the structures and the corresponding operations of linked lists used to store particle information are schematically shown in Fig. 8.9. As seen, each CPU owns one linked list to store information of the particles that are at the corresponding locations to the gas phase. In each linked list, one node represents one particle and consists of two parts: the data part that stores particle information, and the pointer that points to the next node of the particle for traversing. The pointer of the last node always points to the "NULL", implying that the linked list has ended.

Additionally, there are four basic operations for the linked list to adjust the particle's storage location: allocate, free, delete and insert, as depicted in Fig. 8.9. The "allocate" operation is used to allocate new memory to store information about a "new" particle, the "free" operation to free the memory that stores information about an "old" particle, the "delete" operation to disconnect one particle from the linked list, and the "insert" operation to connect one particle at the end of the linked list. When the particle phase is staggered with the gas phase because of relative motion, the information of the particles, whose locations exceed the corresponding location ranges assigned to the present CPUs, will be transferred to and stored in the CPUs with the correct particle location ranges. For example, for CPU A in Fig. 8.9, where one particle is removed, the following sequence of operations is needed:


**Fig. 8.9** Linked lists and operation sequences. Courtesy of Z. J. Zhang [9]

Meanwhile, for CPU B in Fig. 8.9, where a new particle is received, the corresponding operation sequence is:


With this dynamic data structure and the corresponding operation sequences, the information about Lagrangian particles is always stored in the CPUs of the correct Eulerian coordinates, and therefore excessive communication between CPUs when calculating gas-particle interactions is avoided, as shown in Fig. 8.10. Moreover, with the limitation of the global time step by the CFL condition, only "one" particle at most will cross the CPU boundary at every iteration; that is, the information of "one" particle at most will be transferred to the other CPU at one iteration step by the above operation sequence. As a result, the communication cost of the MPI parallel for the gas-particle interaction calculation will be reduced from O(*N*) to O(1) when using this data structure. Here, *N* is the number of particles stored in each CPU.

To demonstrate the MPI parallelization performance with the use of the above linked lists and the corresponding operation sequences, Fig. 8.11 depicts the speedup parameters of a 2D Al-air detonation propagation problem using the CESE method under an Eulerian–Lagrangian framework, which involves approximately 24 million Eulerian meshes and about 75 million Lagrangian particles in the computational domain. This test case is simulated on the Tianhe-2 supercomputer from China with core numbers of 1, 2, 3, 4, 6, 12, 24, 48, 96, 192 and 384. The use of one core means

**Fig. 8.10** Dynamic data structure using the linked list in the Eulerian–Lagrangian framework. Courtesy of Z. J. Zhang [9]

that the simulation is done serially. As depicted in Fig. 8.11, when 384 cores are used, the code using linked lists still has a reasonable parallel efficiency of about 50% for the tested problem, while MPI parallelization under the Eulerian–Lagrangian framework using multi-dimensional static arrays is impossible, even with only 2 cores, as the communication cost is shown to be unacceptably large. This means that the linked list dynamic data structure works well in the MPI implement when solving two-phase detonations under the Eulerian–Lagrangian framework.

To test the performances of CESE method in simulations of two-phase detonation under the Eulerian–Lagrangian framework, Shen et al. [23] simulated the gas–liquid two-phase detonations in the C10H22–O2/air systems with different fuel droplet sizes and equivalence ratios. A deficit in the detonation speed compared to the corresponding purely gaseous one was observed in the gas–liquid suspensions with lean fuel and larger droplet sizes, while an increase in the detonation speed was observed

[9]

**Fig. 8.12** Distribution of particle size near the detonation front in gas–liquid two-phase detonations. Courtesy of H. Shen [23]

with very rich fuel. This special two-phase detonation feature in rich fuel is due to the slow combustion rate of the droplet, leading to a large amount of fuel burnt behind the C–J point, as shown in Fig. 8.12. As a result, the efficient equivalence ratio is reduced and detonation speed shifts to the equivalence ratio direction. The accuracy and robustness of the CESE method applied to two-phase detonation simulations under an Eulerian–Lagrangian framework had been demonstrated.

Notably, simulations of the monodisperse two-phase detonation problems (with single particle size) under the Eulerian–Lagrangian framework always present the same results as those performed under the Eulerian–Eulerian framework. The raising of modelling two-phase detonations under the Eulerian–Lagrangian framework is in fact to simulate the realistic polydisperse two-phase detonations, where a particle size distribution is involved and the particle diameters are in a relatively wide range. One of the frequently used particle size distribution models to describe the polydisperse suspension is the log-normal distribution function, expressed by the number frequency distribution function as follows [24]:

$$f\_n(d\_p) = \frac{1}{\sqrt{2\pi}\sigma\_0} \exp\left[ -\frac{1}{2} \left( \frac{\ln d\_p - \ln d\_{nM}}{\sigma\_0} \right)^2 \right] \frac{1}{d\_p},\tag{8.28}$$

where *dnM* and σ0 are the number median diameter and standard deviation of the distribution, respectively. In practice, the specific particle size distribution of a polydisperse suspension is always given by other two particle size parameters that can be easily measured, namely the volume-average diameter (*d*) and the mass-weightedaverage particle diameter (*dm*). For example, the measured polydisperse parameters of the Al powder tested in the experiment of Zhang et al. [20] are *d* = 2 µm and *dm*  = 3.3 µm. The relationships between (*dnM* , σ0) and (*d*, *dm*) can be obtained through the integrations of Eq. (8.28), as follows

$$\begin{cases} \overline{d} = \left[ \int\_0^\infty s^3 f\_n(s) ds \right]^{1/3} = d\_{nM} e^{\frac{1}{2} \sigma\_0^2} \\ \overline{d}\_m = \left[ \int\_0^\infty s^4 f\_n(s) ds \right] / \left[ \int\_0^\infty s^3 f\_n(s) ds \right] = d\_{nM} e^{\frac{\gamma}{2} \sigma\_0^2} \end{cases} \quad (8.29)$$

Consequently, the corresponding parameters of the above-mentioned Al powder are *dnM* = 1.37 µm and σ0 = 0.5. However, it is a good way to use the parameter set of (*d*, σ0) to discuss the polydisperse suspensions in numerical simulations, because the effects of particle size distribution on polydisperse detonations can be easily discussed by varying σ0. Some specific particle size distributions with fixed *d* = 2 µm and different σ0 values are depicted in Fig. 8.13.

In the followings, the typical Al-Air detonation problems with different particle size distributions, taking as an example, are simulated by the above numerical algorithm under the Eulerian–Lagrangian framework, to show the capacities of the CESE method in solving polydisperse two-phase detonation and show the differences between monodisperse and polydisperse two-phase detonations. The single-step global chemical reaction occurs in the Al-air suspension is

$$2\text{Al}(\text{s}) + 3/2\text{O}\_2(\text{g}) \rightarrow \text{Al}\_2\text{O}\_3(\text{s})\tag{8.30}$$

**Fig. 8.13** Log-normal particle size distributions with *d* = 2 µm. Courtesy of Z. J. Zhang [25]

#### 8.2 Two-Phase Detonations 117

To model the combustion rate of every single Al particle, the surface-kineticoxidation and diffusion hybrid combustion model (originally proposed by Zhang et al. [26]) can be employed. The Lagrangian form of the combustion rate of the *k*th Al particle is given as

$$J\_{pk, \text{Al}} = \pi d\_{pk, \text{Al}}^2 \mathcal{C}\_{\text{O}\_2} \frac{\nu\_{\text{Al}} W\_{\text{Al}}}{\nu\_{\text{O}\_2} W\_{\text{O}\_2}} \cdot \frac{k\_{dk} k\_{\text{x}k}}{k\_{dk} + k\_{\text{x}k}}.\tag{8.31}$$

where the reaction rates for the diffusion-controlled and kinetic-controlled combustion regime are expressed by

$$\begin{cases} k\_{dk} = \frac{\nu\_{\rm O\_2} W\_{\rm O\_2}}{\nu\_{\rm Al} W\_{\rm Al}} \frac{\rho\_{\rm Al} d\_{pk,\rm Al}}{2C\_{\rm total} K d\_{pk0,\rm Al}^2} \left( 1 + 0.276 \text{Re}\_{pk}^{1/2} \overset{1/3}{\text{Pr}} \right) . \end{cases} \tag{8.32}$$

In Eqs. (8.31) and (8.32), νAl and νO2 are the stoichiometric coefficients for Al and O2, respectively; *C*O2 and *C*total are the mole concentrations of O2 and gas mixture, respectively; *Tsk* is the particle surface temperature; and *K*, *k*0 and *Ea* are model constants.

As depicted in Fig. 8.14. Some special two-phase detonation features corresponding to multi-phase interaction can be clearly identified in the detonation front structures of monodisperse Al-air detonation. The first feature is known as "double peaks" in the gas pressure, density and velocity profiles (Fig. 8.14a, c), which is distinctly different from the single peak feature observed in gaseous detonations. It is demonstrated, from the one-dimensional flow theory in gas dynamics that the second peak in the detonation front structures of monodisperse Al suspension is caused by the dominant stage of heat transfer due to intense phase transition (Al evaporation) at a specific location after the shock front. The second feature, shown in Fig. 8.14b, is the plateau of particle temperature due to Al evaporation, which is equal to 2750 K and results in the observed "kink" in the gas temperature profile due to the intense heat transfer between the gas and particles. The third feature, shown in Fig. 8.14c, is the particle velocity lag in the velocity relaxation process, resulting in the alternative forces acting on particles and momentum transfers between the gas and particles. All these two-phase features for monodisperse Al-air detonation are similar with those obtained by Zhang et al. [26] and Teng and Jiang [27, 28].

However, as for the polydisperse detonation with a log-normal particle size distribution of σ0 = 0.5, most features of monodisperse Al detonation, including the double peaks in gas pressure, density, velocity profiles and the kink in gas temperature profile, disappear, and only single peaks of gas quantities exist in the detonation front, which are quite similar to the wave front structures in gaseous detonations. The double peaks in the detonation front is demonstrated theoretically to be attributed to the space-dispersed phase transition processes of particles of different sizes result in an overall moderate heat transfer intensity, which hinders the formation of the heattransfer-dominant stage. These differences between monodisperse and polydisperse

**Fig. 8.14** Comparison of front structures in gas phase of polydisperse and monodisperse Al-air detonations: **a** pressure, **b** temperature, **c** velocity, **d** density, **e** mass fraction of O2. Courtesy of Z. J. Zhang [25]

two-phase detonations are relative to the multiple timescales and length scales in polydisperse suspension with a continuous particle size distribution.

As for the multi-dimensional two-phase detonation features, Fig. 8.15 shows the comparisons of cellular detonation flow fields between monodisperse and polydisperse detonations. As indicated in Fig. 8.15a, typical cellular detonation structures, including pairs of triple points, Mach stems, incident shock waves and pairs of transverse waves, can be observed in both the monodisperse and polydisperse detonation fronts. Moreover, the gas temperature and the O2 species mass fraction distributions for both monodisperse and polydisperse cases are non-uniform behind the detonation fronts, which are characterized by irregular local (high and low) temperature and O2 concentration zones, respectively, as shown in Fig. 8.15b, c. These irregular distributions of the flow field parameters in monodisperse and polydisperse detonations are both caused by the periodical motions of triple points along the detonation fronts. All these features in two-phase detonations are similar to those observed in gaseous detonations, except for the transverse waves. The transverse waves in both the monodisperse and polydisperse detonation fronts are weak and degenerate fairly fast in the rear flows, which are different from the strong transverse waves observed in purely gaseous detonations. According to Zhang et al. [26], these weak transverse waves can be attributed to the slow diffusion-controlled combustion of the majority of Al particles after their kinetic-inductions and a considerable amount of condensed Al oxide formed in the detonation products without contributing to gas pressure.

**Fig. 8.15** Comparison of flow fields in Al-air detonation fronts. Left: monodisperse with *dp* = 2 µm at *t* = 0.6 ms. Right: polydisperse with *d* = 2 µm and σ0 = 0.5 at *t* = 1 ms. Courtesy of Z. J. Zhang [25]

Further, larger detonation cell sizes are expected in polydisperse two-phase detonation compared to monodisperse detonation, since the reaction zone is larger, as shown in Fig. 8.16 by peak pressure contours. The estimated cell sizes of the monodisperse detonation and the polydisperse detonations (σ0 = 0.5 and 0.8) are λmonodisperse = 10.5 ± 0.5 mm, λ0.5 = 13.3 ± 0.8 mm and λ0.8 = 20.0 ± 1.8 mm, respectively. λ0.5 is 27% larger than λmonodisperse, and λ0.8 is even 190% larger. Figure 8.17 plots these detonation cell sizes as a function of the square of standard deviation σ<sup>2</sup> 0 by a logarithmic scatter diagram, together with those of other three polydisperse detonations with different σ0. An approximately linear relationship between ln λ and σ<sup>2</sup> 0 can be captured. Then, linear fitting is employed, and the fitting function appears to be ln λ = 0.9367σ<sup>2</sup> <sup>0</sup>+ 0.0655, implying that the detonation cell size λ of polydisperse suspensions is an exponent function of the square of standard deviation of the distribution σ<sup>2</sup> 0 . This is an important quantitative relation of polydisperse two-phase

**Fig. 8.16** Comparison of cellular detonations (peak pressure contours) of Al-air mixtures: **a**  monodisperse with *dp* = 2 µm, **b** polydisperse with *d* = 2 µm and σ0 = 0.5, and **c** polydisperse with *d* = 2 µm and σ0 = 0.8. Courtesy of Z. J. Zhang [25]

detonation to evaluate the effect of particle size distribution on detonation cell size. From the above results, the capacities of the CESE method in polydisperse two-phase detonation simulations under the Eulerian–Lagrangian framework is demonstrated as well.

#### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **Chapter 9 Other Applications**

The CESE method has been applied to a wide range of scientific and engineering problems since its inception in the 1990s. Although solving CFD problem is the primary goal of the CESE method, this general approach is actually applicable to a variety of PDE systems with physical backgrounds different from fluid dynamics. This chapter mainly introduces the application highlights of CESE in several representative fields.

#### **9.1 Hypersonic Aerodynamics**

Typical high-speed aerodynamic problems can be solved by using the CESE methods. Successful examples of such applications include the supersonic flow over a blunted flat plate [1], Mach reflection of shock waves [2], unsteady viscous flows in rocket nozzles [3], counterflow jets for the reduction of aerothermal load on spacecrafts [4, 5], and the flow over roughness elements in a supersonic boundary layer [6].

With the increasing number of Mars exploration missions in recent years, the research interest of hypersonic aerothermodynamics has been revived. In a typical atmospheric-entry flight, the aircraft will encounter hypersonic gas flow, and the flow field is characterized by strong shock waves and thermochemical non-equilibrium effects. CFD has become a useful tool for in-depth understanding of complex flow physics and for complementary experimental prediction of aerodynamic characteristics.

Using a 2D axisymmetric CESE solver based on hybrid meshes, Wen et al. [7– 9] conducted a systematic study of hypersonic chemically reacting non-equilibrium flows over spheres. Three different working gases: nitrogen, air, and carbon dioxide were used in their simulations. The thermochemical non-equilibrium effects considered therein include the vibrational energy relaxation of gas molecules as well as the dissociation and recombination reactions in the gas flows. The physically consistent coupled vibration–chemistry–vibration (CVCV) [10] two-temperature model was employed. Here, the performance of the CESE method are shown by three different flow cases. The working gas, sphere radius, and the CESE simulation results [9] of the dimensionless shock stand-off distances (defined as the ratio of the stand-off distance to the sphere radius) of each case are listed in Table 9.1, along with the experimental and theoretical results of Wen and Hornung [11]. The CESE results are shown in good agreement with experimental and theoretical results. Further investigation of the shape of shock waves is presented in Fig. 9.1, in which the numerical density contours are overlaid on the experimental finite fringe differential interferograms. The shapes of the simulated bow shocks with the CESE method match the experimental results well.

Another important thermochemical non-equilibrium process in an atmosphericentry flow is the ionization in a high-temperature shock layer. Massimi et al. [12] conducted a numerical investigation of hypersonic ionized air flows over rounded nose geometries, using a 2D axisymmetric CESE solver based on hybrid meshes. The multi-species NS equations were solved with the two-temperature seven-species Park model for ionization reactions. The weakly ionized air flow in the Radio Attenuation Measurement (RAM-C II) flight test at 71 km altitude was simulated and compared with experiments and numerical results of Candler and MacCormack [13]. The CESE results [12] of the maximum electron number densities along the direction normal to the vehicle's surface agree well with the flight test measurements and the numerical results of Candler and MacCormack [13].

**Table 9.1** The working gases, sphere radii, and dimensionless shock stand-off distances (last three columns) in different flow cases


**Fig. 9.1** Comparison of experimental shapes of the bow shocks with the numerical predictions by CESE. Courtesy of C. Y. Wen [9]

#### **9.2 Aeroacoustics**

Thanks to the rapid development of CFD, it is now expected to solve aeroacoustics problems through high-resolution numerical simulations of unsteady compressible Euler or NS equations. However, in many computational aeroacoustic problems, the coexistence of shock waves and small disturbances (acoustic waves) imposes stringent requirements on numerical methods.

The sound–shock interaction problem was used by Wang et al. [14] to examine the accuracy of the CESE method for aeroacoustic problems involving shock waves. Loh et al. [15] further considered three selected problems, namely a linear benchmark problem, instability waves on a free shear layer, and shock–vortex interactions, to illustrate the feasibility of using CESE as a tool for computational aeroacoustics (CAA). Good numerical results with high resolution and low dispersion were demonstrated by their second-order CESE scheme. Additionally, the authors also pointed out that the non-reflecting boundary condition, which plays an important role in CAA, is much simpler to implement in the CESE method than in traditional methods (see [15, 16]). To deal with practical multi-dimensional CAA problems where large disparity in grid cell size is inevitable, Yen et al. [17, 18] modified the Courant number insensitive (CNI)-CESE scheme [17] and the local time-stepping CESE scheme [18] to improve the numerical efficiency. A 3D simulation was performed for the propagation of a Gaussian acoustic pulse and a vorticity wave embedded in a Mach 0.5 mean flow. Once again, the CESE results were in good agreement with the corresponding analytical solution in preserving both the form and amplitude of the waves [17]. Yen [18] further conducted a series of 3D CESE simulations for the Helmholtz resonator problem. By using a second-order CESE scheme, excellent agreement between the linear acoustic theory and the CESE solution was achieved.

The CESE method has also been employed as a numerical tool for CAA problems in practical engineering applications related to engines and aircrafts. To shed light upon an aeroacoustic resonance phenomenon often encountered by convergent–divergent nozzles under transonic conditions, Loh and Zaman [19, 20] developed an axisymmetric NS solver using the CESE method and conducted a numerical investigation. A 3D CESE NS solver was also implemented by Loh et al. [15] and applied to compute the screech noise generated by an under-expanded supersonic jet. Kim et al. [21] simulated the supersonic unsteady flow over an open cavity to investigate the mixing-enhancement and flame-holding capability in the scramjet engine. The CESE method solving 2D NS equations successfully captured the selfsustained oscillations in the supersonic cavity flows. The computed frequencies and amplitudes of the pressure oscillations compare favourably with the theoretical and experimental data. Cheng et al. [22] also studied subsonic and supersonic flows over open-cavity geometries by an unsteady NS solver based on the CESE method. When compared with experimental and analytical data, the generation and propagation of flow-induced acoustic waves were faithfully captured, and satisfactory results were obtained for the oscillation frequency of the dominant mode.

#### **9.3 Solid Dynamics**

In recent years, high-resolution discontinuity-capturing numerical methods originally devised for CFD problems have been utilized to solve nonlinear solid-dynamic problems with Eulerian formulations. Among these methods, the CESE method has attracted considerable attention from researchers in computational solid dynamics because of its simple logic, high accuracy, and ability to capture the behaviours of nonlinear waves.

Wang et al. [23] investigated numerically two high-velocity impact problems in elastic–plastic materials: (1) Taylor copper-bar-impact problem and (2) the penetration of a long rod (tungsten heavy alloy) into a steel target. In this study, the CESE method combined with the level-set technique was adopted to solve Eulerian governing equations that describe the solid dynamics and to trace material interfaces. Chen et al. [24] extended the CESE solver developed by Wang et al. [23] to simulate high-velocity impact problems involving elastic–plastic flows, high strain rates, and spall fractures. Excellent agreement was observed between their simulation of aluminium plates colliding with stainless-steel plates and the experimental data.

As known, the governing equations of elastic waves in solids can be properly written as a set of fully coupled first-order hyperbolic PDEs, including the conservation laws of mass and momentum as well as the rate-type constitutive relations for materials, by treating the density, velocity, and stress components as primitive unknowns. Therefore, the CESE method, which is suitable for hyperbolic systems, can be intuitively applied to simulate linear and nonlinear waves in elastic solids. The CESE simulations of resonant standing waves arising from a time-harmonic external axial load and compression waves arising from a bi-material collinear impact was demonstrated by Yu et al. [25]. Another study on the propagation and reflection of extensional waves in an abruptly stopped elastic rod using the CESE method was carried out by Yang et al. [26]. Moreover, Chen et al. [27] performed CESE simulations to study planar-wave expansion from a point source in an anisotropic solid of cubic symmetry (gallium arsenide) and Yang et al. [28] extended the CESE simulations to stress waves in solids of hexagonal symmetry, illustrating wave propagation in a heterogeneous solid composed of three blocks of beryl with different lattice orientations. Lately, Lowe et al. [29] conducted a comprehensive study of nonlinear longitudinal waves (including both weak and strong shocks, rarefactions, and contact discontinuities) in tapered elastic rods using the CESE method as a numerical tool. In the above studies, the CESE simulations effectively captured waves in a variety of solid materials and provided results consistent with the available analytical solutions.

#### **9.4 Magnetohydrodynamics**

The plasma flows in aerospace applications and astrophysics can be described and solved generally by the magnetohydrodynamic (MHD) equations, which combine the NS equations for fluid dynamics and the Maxwell equations for electromagnetics. Due to the divergence-free constraint usually imposed for the magnetic field *B* in the computational MHD problems, some special treatments are added into the MHD numerical schemes to enforce ∇ · *B* = 0. Since the proposal of the CESE method, researchers have been attempting to solve MHD problems accurately and efficiently with the CESE method.

Zhang et al. [30, 31] used the CESE method to study MHD benchmark problems, including an MHD shock-tube problem and an MHD vortex problem. Their CESE results without additional treatments for ∇ · *B* = 0 compared favourably with previously reported reference solutions. In fact, Zhang et al. [30, 31] performed the simulations with the baseline CESE scheme and the CESE scheme in conjunction with a special treatment to maintain ∇ · *B* = 0. Nevertheless, no obvious difference was observed.

Feng et al. [32] developed a numerical platform based on the CESE method to investigate solar–interplanetary physics and space weather. Three dimensional MHD equations were solved with the adaptive mesh refinement (AMR) to better resolve flow features that have spatial scales many orders of magnitude smaller than the vast size of solar–interplanetary space. This CESE-based MHD-simulation approach was validated through the numerical study of solar corona and solar wind and comparison with observation data. Recently, Yang et al. [33, 34] extended this CESE-based MHD solver to a high-order version [33] and an upwind version [34] based on the CESE modifications of Shen et al. [35] and Shen and Wen [36], respectively. Numerical tests of MHD-vortex and MHD-blast-wave problems demonstrated the accuracy improvement of the CESE results by applying both the high-order and the upwind extension. Furthermore, a new strategy to keep the magnetic field fundamentally divergence-free was proposed by Yang et al. [33], taking the advantage of the CESE method that the spatial derivatives of the magnetic field are treated as marching variables in the algorithm.

Notably, the successful combination of MHD modelling and CESE simulation was elaborated in a recent monograph of Feng [37]. For more details about the CESE method for MHD, readers are recommended to refer to Feng [37] and the references therein.

#### **References**


**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **Chapter 10 Summary**

In this book, the space–time CESE method is introduced as a novel and distinctive numerical approach to solving conservation equations in engineering sciences. Also, from a historical perspective, the development of the CESE method since the 1990s has been reviewed and the remarkable improvements and extensions of the CESE method in the past few years has been summarized. Various applications of the CESE method have been presented to emphasize its numerical performance under different scenarios, including many research topics in engineering, such as compressible multifluid flows and detonations.

The most elegant part of the CESE method might be the non-dissipative *a* scheme described in Chap. 2, where the unknown variable *u* and its spatial derivative *ux*  can be obtained directly from the conservation laws. However, for the purpose of shock capturing, two families of CESE schemes have been developed on the basis of *a* scheme. The first family is called the central CESE scheme, including the *a–*α and CNI schemes. In these schemes, the spatial derivatives are updated by specially designed procedures, where artificial dissipation can be added. The central CESE schemes avoid any characteristic-based techniques (i.e., they are free of Riemann solver). The second family is the upwind CESE schemes, in which the upwind flux solver (e.g., Riemann solver) can be utilized, but it just plays a supporting role. By combining the ideas in the original CESE method and the conventional upwind FVM, the upwind CESE scheme managed to retain the forms of discretized equations for *u* and *ux* in the *a* scheme, while the numerical dissipation can be introduced in a reasonable manner. It is worth noting that the dissipation of the upwind CESE scheme is naturally insensitive to the CFL number.

Two important issues about the CESE method, namely the CESE schemes on unstructured meshes and the high-order CESE schemes, have been discussed in Chaps. 4 and 5, respectively. The CESE method is naturally multi-dimensional and compatible with unstructured meshes. A high-order extension of the CE/SE method in both space and time can be achieved by storing and updating high-order derivatives (e.g., *uxx* and *uxxx*) at each solution point, without enlarging the stencil or adding stages in time integration.

The CESE method has low dissipation and high compactness. When applied to the simulations of complex physical processes, the CESE method can catch shock waves, contact discontinuities, fine structures, and small disturbances, with high resolution and strong robustness. Therefore, the CESE method demonstrates good performances in the numerical simulations of wave-propagation problems (e.g., shock waves, acoustic waves, detonation waves, stress waves in solid, and the electromagnetic waves), interfacial instabilities, and shock–bubble/droplet interactions. In many hot research areas including hypersonic aerodynamics, shock dynamics, detonation, aeroacoustics, MHD, and solid dynamics, the CESE method proves to be suitable and shows a good development prospect.

In recent years, an important topic in CESE research has been its application to highly nonuniform and high-aspect-ratio computational meshes, which are often required in simulations of boundary layers. In particular, the CESE method has been recently applied to aerodynamic heating problems and direct numerical simulation of laminar or turbulent flows. It is also worth noting that the CESE schemes are mainly applied to unsteady problems due to their explicit nature in time with the associated restriction on time step size. In order to solve steady problems efficiently, the extensions to implicit time-stepping CESE schemes are under developing. Further analyses of the numerical properties of CESE schemes will be very useful to solidify their mathematical foundation. The modified equation analysis and the spectral property analysis for the upwind CESE scheme are in progress now.

Along with the continuous improvement of high-performance computing, numerical methods emerge one after another. Many challenging numerical simulations can be accomplished today. However, none of the existing numerical methods can be accurate, robust, and efficient under all situations. In this sense, the CESE method may not be superior to the existing methods, but it does provide an alternative approach which deserves more evaluations and further investigations. The authors wish this book can familiarize the readers with the CESE method.

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

## **Appendix**

1. The *a*-α CESE solver for solving the 1D scalar convection problem.

```
#include <stdio.h> 
#include <math.h> 
#include <stdlib.h> 
#include <time.h> 
#define end_time 2.0 // Computing time. 
#define L1-1.0 // Left end of the domain. 
#define L2 1.0 // Right end of the domain. 
#define N 200 // Mesh number. 
#define Sf 0.8 // CFL number. 
#define MAX(x,y) (((x)>(y))?(x):(y)) 
#define MIN(x,y) (((x)<(y))?(x):(y)) 
double U1[N+3],U2[N+3],Ux1[N+3],Ux2[N+3],Ut[N+3],F[N+3],Ft[N+3]; 
double const dx=(double(L2-L1))/double(N); //Mesh size 
double avg(double x1,double x2) 
{ 
      double x,aa; 
      aa=2.0;//6*fabs(x1-x2)/(fabs(x1)+fabs(x2)+1e-10); 
      x=(pow(fabs(x1),aa)*x2+pow(fabs(x2),aa)*x1)/ 
            (pow(fabs(x1),aa)+pow(fabs(x2),aa)+1e-20); 
      return x; 
}// Compute weighted average of x1 and x2. 
double CFL(double U[N+3]) 
{ 
      int i; 
      double maxvel,vel; 
      maxvel=1e-10; 
      for(i=1;i<=N+1;i++) 
      { 
            vel=1.0; 
            if (vel>maxvel)maxvel=vel; 
      } 
      return Sf*dx/maxvel; 
}// Compute timestep.
```

```
© The Editor(s) (if applicable) and The Author(s) 2023 
C.-Y. Wen et al., Space–Time Conservation Element and Solution Element Method, 
Engineering Applications of Computational Methods 13, 
https://doi.org/10.1007/978-981-99-0876-9
```
133

```
void Bound (double U[N+3],double Ux[N+3]) 
{ 
      U[0]=U[N]; 
      U[N+2]=U[2]; 
      Ux[0]=Ux[N]; 
      Ux[N+2]=Ux[2]; 
}// Apply boundary condition. 
void U2Fu(double U[N+3],double Fu[N+3]) 
{ 
      int i; 
      for(i=0;i<N+3;i++) 
      { 
            Fu[i]=U[i]; 
      } 
}// Compute flux. 
void ComputeTemporalDerivatives(double U[N+3],double Ux[N+3]) 
{ 
      int i; 
      for(i=0;i<=N+2;i++) 
      { 
            Ut[i]=-Ux[i]; 
            Ft[i]=Ut[i]; 
      } 
}// Compute the temporal derivatives of U and F. 
void CESE_TimeMarch(double U[N+3],double Ux[N+3],double 
Fu[N+3],double 
Fut[N+3],double Unew[N+3],double Uxnew[N+3],double dt,int ishalf) 
{ 
      int i,I; 
      double Uxplus,Uxminus,Uleft,Uright,Fleft,Fright; 
      for(i=1-ishalf;i<=N+1;i++) 
      { 
            I=i+ishalf; 
            Uleft=U[I-1]+dx*Ux[I-1]/4; 
            Uright=U[I]-dx*Ux[I]/4; 
            Fleft=Fu[I-1]+dt*Fut[I-1]/4; 
            Fright=Fu[I]+dt*Fut[I]/4; 
            Unew[i]=0.5*(Uleft+Uright+dt/dx*(Fleft-Fright)); 
            Uxminus=2*(Unew[i]-U[I-1]-dt*Ut[I-1]/2)/dx; 
            Uxplus=2*(U[I]+dt*Ut[I]/2-Unew[i])/dx; 
            Uxnew[i]=avg(Uxminus,Uxplus); 
      } 
}// Scheme marching from integer-step to half-step, or from half-
step to integer-step. 
void Initialize(double U[N+3],double Ux[N+3]) 
{ 
      int i; 
      for(i=0;i<=N+2;i++) 
      {
```

```
if(i > =(N/4+1)&&i<=(3*N/4+1)) 
            { 
                   U[i]=1; 
            } 
            else 
            { 
                   U[i]=0; 
            } 
            Ux[i]=0; 
      } 
}// Initialization. 
void Output(double U[N+3],double t) 
{ 
      int i; 
      double x; 
      FILE *fp; 
      fp=fopen("result.txt","w+"); 
      for(i=1;i<=N+1;i++) 
      { 
            x=L1+(i-1)*dx; 
            fprintf(fp,"%20.5e\t%20.5e\n",x,U[i]); 
      } 
      fclose(fp); 
}// Output result. 
void CESE_Solver() 
{ 
      int Nstep=0; 
      double dt,t=0; 
      Initialize(U1,Ux1); 
      while(t<end_time &&Nstep>=0) 
      { 
            dt=CFL(U1); 
            t=t+dt; 
            Nstep++; 
            printf("Nstep=%d,dt=%e,t=%e\n",Nstep,dt,t); 
            U2Fu(U1,F); 
            ComputeTemporalDerivatives(U1,Ux1); 
          CESE_TimeMarch(U1,Ux1,F,Ft,U2,Ux2,dt,1);//the 1st half 
time step 
            U2Fu(U2,F); 
            ComputeTemporalDerivatives(U2,Ux2); 
          CESE_TimeMarch(U2,Ux2,F,Ft,U1,Ux1,dt,0);//the 2nd half 
time step 
            Bound(U1,Ux1); 
      } 
      Output(U1,t); 
}// a-α CESE solver. 
int main(void) 
{ 
      CESE_Solver(); 
      return 1; 
}
```

```
2. The a-α CESE solver for solving the Sod-shock problem.
```

```
#include <stdio.h> 
#include <math.h> 
#include <stdlib.h> 
#include <time.h> 
#define end_time 0.1 // Computing time. 
#define L1 0.0 // Left end of the domain. 
#define L2 2.0 // Right end of the domain. 
#define L0 1.0 
#define N 200 // Mesh number. 
#define Sf 0.8 // CFL number. 
#define GAMA 1.4 // Ratio of specific heats. 
#define rou1 1.0 
#define u1 0.0 
#define p1 1.0 
#define rou2 0.125 
#define u2 0.0 
#define p2 0.1 // Initial conditions 
#define MAX(x,y) (((x) > (y))?(x):(y)) 
#define MIN(x,y) (((x)<(y))?(x):(y)) 
double U1[N+1][3],U2[N+1][3],Ut[N+1][3],F[N+1][3],Ux1[N+1][3], 
Ux2[N+1][3], Ft[N+1][3]; 
double const dx=(double(L2-L1))/double(N); // Mesh size. 
double CFL(double U[N+1][3]) 
{ 
     int i; 
     double maxvel=1e-4,vel,p,u; 
     for(i=1;i<=N-1;i++) 
     { 
           u=U[i][1]/U[i][0]; 
           p=(GAMA-1)*(U[i][2]-0.5*U[i][0]*u*u); 
           vel=sqrt(GAMA*p/U[i][0])+fabs(u); 
           if (vel > maxvel)maxvel=vel; 
     } 
     return Sf*dx/maxvel; 
}// Compute timestep. 
double avg(double x1,double x2) 
{ 
     double x,aa=1.0; 
     x=(pow(fabs(x1),aa)*x2+pow(fabs(x2),aa)*x1)/ 
            (pow(fabs(x1),aa)+pow(fabs(×2),aa)+1e-20); 
     return x; 
}// Compute weighted average of x1 and x2. 
void Init() 
{ 
     int i,k; 
     double x; 
     for(i=0;i<=N;i++) 
     { 
           x=L1+i*dx; 
           if(x<=L0)
```

```
{ 
                   U1[I][0]=rou1; 
                   U1[I][1]=rou1*u1; 
                   U1[I][2]=p1/(GAMA-1)+rou1*u1*u1/2; 
            } 
            else 
            { 
                   U1[I][0]=rou2; 
                   U1[I][1]=rou2*u2; 
                   U1[I][2]=p2/(GAMA-1)+rou2*u2*u2/2; 
            } 
      } 
      for(i=0;i<=N;i++) 
      { 
            for(k=0;k<3;k++) 
            { 
                   Ux1[i][k]=0; 
            } 
      } 
}// Initialization. 
void Bound (double U[N+1][3]) 
{ 
      for(int k=0;k<3;k++) 
      { 
            if(k==1) 
            { 
                   U[N][k]=0; 
                   Ux1[N][k]=Ux1[N-1][k]; 
                   U[0][k]=0; 
                   Ux1[0][k]=Ux1[1][k]; 
            } 
            else 
            { 
                   U[N][k]=U[N-1][k]; 
                   Ux1[N][k]=0; 
                   U[0][k]=U[1][k]; 
                   Ux1[0][k]=0; 
            } 
      } 
}// Apply boundary condition. 
void U2F(double U[N+1][3],double F[N+1][3]) 
{ 
      int i; 
      double u,p; 
      for(i=0;i<N+1;i++) 
      { 
            u=U[i][1]/U[i][0]; 
            p=(GAMA-1)*(U[i][2]-0.5*U[i][1]*U[i][1]/U[i][0]); 
            F[i][0]=U[i][1]; 
            F[i][1]=U[i][0]*u*u+p; 
            F[i][2]=(U[i][2]+p)*u; 
      } 
}// Compute flux.
```

```
void ComputeDerivative(double U[N+1][3],double Ux[N+1][3]) 
{ 
      int i; 
      double u; 
      for(i=0;i<=N;i++) 
      { 
          u=U[i][1]/U[i][0]; 
          Ut[i][0]=-Ux[i][1]; 
          Ut[i][1]=-(GAMA-3)*u*u*Ux[i][0]/2-(3-GAMA)*u*Ux[i][1]
           Ut[i][2]=-((GAMA-1)*u*u*u-GAMA*u*U[i][2]/U[i][0])* 
Ux[i][0]-(GAMA*U[i][2]/U[i][0]-3*(GAMA-1)/2*u*u)*Ux[i][1]-
GAMA*u*Ux[i][2]; 
           Ft[i][0]=Ut[i][1]; 
           Ft[i][1]=(GAMA-3)*u*u*Ut[i][0]/2+(3-GAMA)*u*Ut[i][1] 
+(GAMA-1)*Ut[i][2]; 
           Ft[i][2]=((GAMA-1)*u*u*u-GAMA*u*U[i][2]/U[i][0])* 
Ut[i][0]+(GAMA*U[i][2]/U[i][0]-3*(GAMA-1)/2*u*u)*Ut[i][1] 
+GAMA*u*Ut[i][2]; 
      } 
}// Compute the derivatives of F. 
void CESE_TimeMarching(double U[N+1][3], double Ux[N+1][3], 
double Fu[N+1][3],double Fut[N+1][3], double Unew[N+1][3], double 
Uxnew[N+1][3], double dt, int ishalf) 
{ 
      int i,j,I; 
      double Uleft,Uright,Fleft,Fright,Ux_minus,Ux_plus; 
      for(i=1-ishalf;i<N;i++) 
      { 
            I=i+ishalf; 
            for(j=0;j<3;j++) 
            { 
                  Uleft=U[I-1][j]+dx*Ux[I-1][j]/4; 
                  Uright=U[I][j]-dx*Ux[I][j]/4; 
                  Fleft=F[I-1][j]+dt*Ft[I-1][j]/4; 
                  Fright=F[I][j]+dt*Ft[I][j]/4; 
                  Unew[i][j]=0.5*(Uleft+Uright+dt/dx*(Fleft
                Ux_minus=2*(Unew[i][j]-U[I-1][j]-dt* 
Ut[I-1][j]/2)/dx; 
               Ux_plus=2*(U[I][j]+dt*Ut[I][j]/2-Unew[i][j])/dx; 
                Uxnew[i][j]=avg(Ux_minus,Ux_plus); 
            } 
      } 
}// Scheme marching from integer-step to half-step, or from half-
step to integer-step. 
void CESE_1D_Solver() 
{ 
      int i,comput_num; 
      double dt,t=0; 
      FILE *fp; 
      comput_num=0; 
      Init();
```

```
while(t<end_time) 
      { 
            comput_num++; 
            dt=CFL(U1); 
            t=t+dt; 
            U2f(U1,f); 
            ComputeDerivative(U1,Ux1); 
            CESE_TimeMarching(U1,Ux1,F,Ft,U2,Ux2,dt,1); 
            U2f(U2,f); 
            ComputeDerivative(U2,Ux2); 
            CESE_TimeMarching(U2,Ux2,F,Ft,U1,Ux1,dt,0); 
            Bound(U1); 
      } 
}// a-α CESE solver. 
void Output(double U[N+1][3]) 
{ 
      int i; 
      double x; 
      FILE *fp; 
      double rou,u,p; 
      fp=fopen("result.txt","w+"); 
      for(i=0;i<N+1;i++) 
      { 
            x=L1+i*dx; 
            rou=U[i][0]; 
            u=U[i][1]/U[i][0]; 
            p=(GAMA-1)*(U[i][2]-0.5*U[i][0]*u*u); 
            fprintf(fp,"%20.5e\t%20.5e\t%20.5e\t%20.5e\n", 
x,rou,u,p); 
      } 
      fclose(fp); 
}// Output result. 
int main() 
{ 
      CESE_1D_Solver(); 
      Output(U1); 
      return 1; 
}
```